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FOREWORD 


This  computer  program  User's  Manual  was  prepared  by 
the  Los  Angeles  Division  of  Science  Applications, 
Incorporated,  for  the  Vehicle  Dynamics  Division  of  the 
Air  Force  Flight  Dynamics  Laboratory,  Wright-Patterson 
Air  Force  Base,  Ohio.  The  computer  programs  were 
developed  under  Project  1370,  "Dynamic  Problems  in 
Flight  Vehicles",  Task  137004,  "Design  Analysis", 

Contract  f 3 36 1 5— 74 — C— 30 9 4 . James  J.  Olsen  and  later 
Lt.  William  L.  Holman  (AFFDL/FYS)  were  the  Air  Force 
Task  Engineers. 

R.  M.  Traci  was  the  principal  investigator  for 
the  study  and  J.  L.  Farr,  Jr.,  developed  the  computer 
programs  described  in  this  report.  Consultant  E.  D.  Albano 
contributed  to  the  development  and  implementation  of 
the  numerical  method. 

The  authors  submitted  this  report  in  July  1975 
for  publication  as  an  AFFDL  technical  report. 

Other  reports  prepared  and  submitted  under  the 
aforementioned  contract  are:  AFFDL-TR-74-37,  "Small 

Disturbance  Transonic  Flows  about  Oscillating  Airfoils," 
AFFDL-TR-74-135,  "Computer  Programs  for  Calculating 
Small  Disturbance  Transonic  Flows  about  Oscillating 

AFFDL-TR- 75-100 , "Small  Disturbance  Transonic 
Flows  about  Oscillating  Airfoils  and  Planar  Wings." 
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1-°  INTRODUCTION 


disturbanoeep0?entiarf™wT?he«v  fDUTR?  lmPlem6"t  a small 
steady  transonic  flow  about'thii^of0*  three'di'"ensioual  un- 
harmonac  oscillation.  The  theorv^.T  "i"95  und«going 

flnlnear  system  “an  be  obtained  bv  e ba?®d  ?"  the  fact  that 
flow  as  a small  perturbation  to  the  ™Sld<Tin9  the  “"steady 
I he  perturbation  exoaniinn  non-uniform  mean  flow 

developed  with  different  emn£Pf°aC?  h?s  recently  been 
hy  the  present  authors  1 ' ^and^hlS  i4dependent  studies 
tions  of  the  theory  and  nuneSLf  l Detailed  descrip- 

the  three-dimensional  version  S?  thfn  method  used  in 
Vsers  roanual,  are  presento??  Programs,  documented 
method  is  a generalize*-  i ln  References  1-3  Th#» 

dimensional  airfoils  in  Refllenc* S de®Sribed  f°r  two- 
programs  documented  in  Reference  21  *£!?  °? . fche  2-D  computer 
the  present  phase  of  rew.rJS  5 2*  . The  final  report  of 
and  presents  some  illustrative  resSlts?8  th*  aeneralization 

per turbationPpotential°f unction^s  appro*ch  Vsed'  the 

increasing  powers  of  e QTT1  i 1 s exPanded  in  a series  of 
of  the  ampl?tule  o!  an  ™*a^K?rK,,hleh  ls  a 
he  resulting  expansion  of  the^um**-**1^31106  t0  the  boundary, 
results  in  a sequence  of  partial  Hiffdy  Potential  equation 
the  perturbation  potentials.  The  sfrofh tlal  “<I“ations  for 
the  usual  nonlinear  steady  tran^oni „r0t?  order  equation  is 

mixed  elliptic/hyperbolic^vpe  and  ? Potential  equation  of 

cSL"1Xe?hfi^erencln9-  "olaxatiSn  procedure  tt  1°^™  USi"7 * 9 
v-oie  . The  first  order  nrn^.,  procedure  of  Murman  and 

and  for  harmonic  boundary  disturbance^*1  ?quation  is  linear 
elliptic/hyperbolic  type?  de^ndinV  18J!lso  of  the  mixed 

Ifc  iS  s°lved  in  TDUTRN  usina  at  uP°n  the  steady  solution, 
used  in  TDSTRN.  9 th  same  numerical  technique  as 

tion  are  discussed  in^he*  following  g0mpVter  Program  opera- 
turbation  theory  and  numerical " tions‘  The  sniall  _ 

summarized  in  Sections  2 0 and  !°iutlon  Procedure  are 
tion  of  the  Program- V ]og?caJ  op;^,rSPeC^iVely-  A ^scrip- 
description  are  given  inSeStiS  4 n1  eand.a  brief  subroutine 
complete  description  of  the  prooram°i  Sfctl°n  5*°  presents  a 
values  for  various  control  vLriab?* a l pU^'  uith  su9gested 
Section  6.0  describes  the  Drnor*m  ' and  Program  output, 
tions  for  making  effective  u«ip  nf  ^a^e  and  includes  sugges- 
which  exercise  all  prooram  tbe  Programs.  Sample  cases 

7-0  With  a completePspeclficatlonSofr!i1?r?8ented  in  9"“«°n 

Finally,  complete  FORTMrj  fwJ  input  and  sample 

TDUTRN  are  presented  in  the  appellees  ?S  °f  TDSTRN  and 


i 


• • -u.  TOI1  ^ 

tool  for  all  speed^ngls^nrL^b^0  Principal  analytical 
tant  in  the  transonic  speed  rannf  Pecome  increasingly  impor- 

used  ^^L^ls9  ^ J*“  SS- 

an“edriCal  • stLSriSe°ti0n  ^e  required 

are  described  in  Section  3 n ;?e.  steady  and  unsteady  system* 

TDSTRNPr°^ramniin9  practice»  the  FORTRAN ted  foll°wing 

TDSTRN  and  TDUTRN  are  descriptive  of^LVa^a?leS  Used  in 

Ptive  of  the  physical  variables. 


M«'U 


F,GURE  L SCHEMATIC  of  T«EE-D,HEH8.«,AL  Planar  WlNG 

(two-diIensionairorfJa^reSt  is  the  flow  about  an  airfoil 
lating  with  various  flexible''^9  ( three-dimensional ) oscil- 

fieltdescthranS°niC  iS'&Eft  de9rees  ^ freedom 

ield  schematic  and  coordinate  *ff° 11  geometry,  flow- 

foil 'T6  ReCtangular  coordinates  f£lt10?  are  given  in  Figure 
foil  leading  edge  with  oriqW  fiX'yfz)  are  fixed  to  the  air- 
are  the  frees tream  velocit?  *lng  root  and  u,  M , a 

spectively.  The  airfoU  has  ”a?*n?mber  and  sound  speeSre- 

attfek  rtmaXimUm  thfcknesshdivided1b^ni?*  r!!tio  6'  wh^ch  is  the 
that  6<<iand^a  semi“sPan  b (ZSPAN)  The  ^hord  Cf  and  angle  of 
' 6<<1  and  ‘ iS  °f  sme  order  o^Kst^ 


i 


Also,  the  oscillatory  motion  of  the  airfoil  is  assumed  to  be 
described  by  a small  non-dimensional  displacement  and 

a reduced  frequency  k=^oc/U  based  on  airfoil  chord  where  w is 
the  frequency  of  oscillation. 

Assuming  inviscid,  isentropic  flow,  the  problem  can 
be  reduced  to  the  solution  of  a single  equation  for  a velo- 
city potential  plus  the  tangency  boundary  condition  on  the 
airfoil  surface.  As  is  well  known,  the  derivation  of  a small 
disturbance  theory  for  transonic  flows  requires  a singular 
perturbation  approach.  The  following  scaling  is  thereby 
introduced: 


x ^ 


t = 


[{1+Y)6mM 


2/3 


^ 1/3 

z = [(l+y)6t£] 


” t 
c 


z 

c 


(1) 


and  the  total  potential  is  expanded  about  the  uniform  flow: 
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Ucx 
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l(l+Y)M’] 
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'Xj  .'Xj  'Xj  'Xj  *\i . 

4>  (x,y,z,t) 


(2) 


Retaining  all  terms  of  leading  order  in  to  total  potential 
equation  and  boundary  conditions  results  in  the  following 
form  for  the  unsteady  small  disturbance  system. 


{ K— = 2$'VVi  + p- 

x xx  ryy  zz  xt  ft  tt 


(3) 


where  the  transonic  similarity  parameters  are: 


(1-M2 ) 


K = 


1 (1+Y)$M2] 


2 1 27T 


ft  = 


K 

[ (1+Y)<5M2]  2/3 
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with  boundary  conditions: 


J: 


% = 
y 


r—  + 2.  _ji » c .'v  % 'v 

n at J fu,« x,z,t) 


on  y = 


0 < x < 1 

- 'b 

0 s z r b 


(4) 


r- 

h + lk 

o 

II 

c>>, 

c 

o 

m 

o 

II 

h + ^ + 

1 

2 

^ 4 0 as  x7+y 

x > 1 

0 5 2 < b 


(5) 


(6) 


(7)rbe^f  on  the  SSS“£,a jrfoil  *»aPe  Unction  (Equation 
where  |]  denotes  a ^mp  " \}Z^JultloeB  «*P«tively,  aSd 
and  0+.  it  is  noted  that  the ?X??ed  quantity  between  y-0" 
tion  (Equation  4)  and  the  Kutta  tangency  boundary  condi- 

aPPUed  in  the  small^^V^T™  ^0“?“°"  5>  are 

the  un  steady ^a i rfoi 1 p?Xe,Hn  ^h  Pr°Vi??S  a formulation  of 
includes  flowfields  with  shocks  Cer?n~U?ear  d?main'  which 
system  are  underlined  as  thee  ml  f-ertain  terms  in  the  above 
quenoy  [kM)(«*/i,j  1 ‘ "ay  be  omitted  for  a low  fre- 

TDUTRN  includes  the  low  ""quench  Pfes!?t  version  of 

^"otal  frequency  formulation  U0p^lf°?  T1*-0!:  <I(''PT=0>  or 
at  the  discretion  of  the  user.  *5PT~1)  and  either  can  be  used 

linear  sJsteT^  a^lve"  (Eaua^"  5°?,  s?lvl"»  *h.  non- 
turbation  potential  function^  terms d"  to  e)tPand  the  per- 
disturbance  c«l.  From  this  point™  alVuid“r?M  ^be^ 


vnr?*h?««ithuthe  understandin9  that  all  variables  are  scaled 
treated?  Harinonic  boundary  disturbances  are  explicitly 


f(x,z,t)  a f (x,z)  + cf  (X,z)eint 
o e 


(7) 


and  the  perturbation  potential  is  expanded  as  follows: 


$ (x,y,z,t) 


4>°(x,y,z)  + e (x,y,z)eint  + ... 


(8) 


Substituting  this  into  the  perturbation  potential  equation  plus 
boundary  conditions  and  combining  terms  results  in  the  follow- 
lng  pair  of  boundary  value  problems  for  <p°  and  <J> 1 respectively. 
(In  the  following  text,  the  superscript  has  been  dropped  from 
<P  • ) 


<k-*x»*xx  + ♦ ♦;«  - ° 

* fg(x,z),  on  y * ± 0 (osxil 

(0  * z * b 

UXJ  “0,  on  y « 0 |x>l 

( 0 < z < b 

(4>^)  +(4»y)  + (<J>°)  ■*  0 as  x2+y2+z2  » 

and 

(K'*xUxx+*yy+*z2  -(^x+2ifl)*x  + iSl  - 0 

♦v  * £^+ikfr,  on  y = i 0 (osxsl 

( 0 i z s b 

jx  > 1 
( 0 $ z £ b 

x2+y2+z2-*'®> 


l$x  + lilt)  “0#  on  y n o 

2 2 2 
(4>x)  + (<J>y)  + (<}>z)  -*  0,  as 


(9) 


(10) 
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System  9 is  recognized  as  the  usual  formulation  for  steady 
transonic  flow  and  system  10  is  the  formulation  for  the  un- 
steady perturbation  thereof.  It  is  noted  that  the  governing 
equation  for  <f>  is  linear  but  of  the  same  mixed  elliptic/ 
hyperbolic  type  as  the  steady  solution.  It  is  also  noted 
that  $ is  in  general  complex  thereby  permitting  phase  shifts 
between  field  quantities  and  the  boundary  disturbance.  As 
before,  underlined  terms  in  system  10  are  neglected  for  a 
consistent  low  frequency  approximation.  Also  for  two  dimen- 
sional airfoil  sections,  the  z dependence  on  all  quantities 

and  the  <j>  terms  in  the  equations  are  neglected, 
z z 

The  main  physical  quantities  of  interest  are  the  pres- 
sure coefficient  and  airfoil  force  coefficients.  The  pressure 
coefficient,  defined  in  the  usual  manner,  is  given  by: 


C 

P 


2/3 


l (1+y)mM 


(CO  + eC  elQt) 

1/3  P P 


(ID 


where  the  steady  and  unsteady  scaled  pressure  coefficients  are 
given  to  leading  order  in  the  small  disturbance  approximation 
by: 


Cp  = "2*x'  Cp  * "2(*x  + ^ 


(12) 


The  formulations  of  the  boundary  value  problems  are 
essentially  complete  with  the  exception  of  the  practical  matter 
of  setting  the  boundary  conditions  away  from  the  airfoil,  which 
depends  on  the  particular  problem;  subsonic  or  supersonic  free 
field,  wind  tunnel  wall  etc.  Asymptotic  far  field  solutions 
to  Equations  10  have  been  developed  for  two-dimensional  sub- 
sonic or  supersonic  free  air  or  wind  tunnel  flows  and  for  three- 
dimensional  subsonic  flow.  These  solutions  are  described  in 
the  present  three-dimensional  subsonic  free  air  version  of  the 
computer  programs. 


^ • o NUMERICAL  solution  method 


prob  JUs  forethea«f»^Uti°S  Prooedure  for  the  boundary  value 

§Hfi 

to  the°governing^i0nS) 

ss^rg?  ssi»H  ^*s- 

r general  lifting  airfoils  developed  by  Krupp6 . 

two-dimens iona^air foil  i*?""*  and  s°l“tion  method  to 
are  interes?^ 

must  be  considered!  flotr4* 

ranges^  K^rtSi'  SpCed=  tha"  in  the  other  speed 

realistic  three-dimensional  c^lluons^prlctJIal^r 
modern  computers  and  it  is  the  purpose  If  ^is  secUoi  ?I 

durertoeperIir"IaIa1Iiu?i:„£stlJIeSSi!?i,  ?°lUti°n  Pr°Ce* 

s:sts“k£™H;^ 
s-sS-O^-v""  » 'ssiaiT 

solutions  to  the  go!Irf  SSjTb^i"^.^^0 

“'=  frce-at-a-  flow  ly  f°or 

oerturh^Ln1  4>k  d by,the  Present  authors  for  the  unsteady 
perturbation.  These  solutions  are  summarized  in  Section  i ? 

SS^?,»s;?cLtSS,4-Bd 

Extensions  of  the°sol"iII6™£hId°fIr'1hI  Ils^dy^Irturba- 
tion  parallel  the  steady  method  and  these  arl  nlw  dl^i^: 
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3.1 


Figure  2.  Schematic  of  Numeral  Solution  Domain 

most  stralghtforwaJdeexte"sionUofrthalt“Cheme  constitutes 
previously  described  in  de’aH  in  two-dimensional  meth 
matically  in  Pigure  2,  a cubic  ree^T"  l-  *<■  shown  s, 
extent  with  uneven  grid  line  sofci^n  ^Ur  m?sh  of  finite 
solution  space.  The  nri^  2 _ ^ °verlayed  on  the  3- 

expanded  oSt  to  the^a^So^^"^^  ne“  the  »ir£n 
difference  equations  are  identical  to  5k  9rid’  The 
dimensional  versions1'1  with  thf  corresponding  twe 

ence  form  for  4 given  by:  addltion  of  a centered  dif 


22  . 

.k 


w,k*4,k-i» 


1 

Az, 


i»j»k+l  ^i,j,k^ 


(13) 


Az <*,  < 


Jk-1 


»jrk  ■ 4>i,j,k-l) 


Figure  3,  Schematic  of  Difference  Scheme 

is  shown  LhcZticnllyainSFi tatovhrteThimenS10nal  8chl-'mc 
elliptic  on  hyperbolic  natu«  of  Ll  \ The  tests  for  the 
centered  difference  form  of  tit-f')*  eqUatl?"  are  Mde  on  the 
value  of  this  coefficient  the  x and  dePending  on  the 

or  backward  differenced  as  in  th!  ? n 8 °f  * are  Ce«tered 
point  operator  is  uSed  iS  to?h  S ? !86*  A Paral*>lic 
and  a shock  point  operator  is  used  f?d  unsteady  schemes 

scheme  for  the  steady  po?en£ia!f  finite  diffe^nce 

for  each  column6 ix^z^^onstantsf6^6”^  equations  are  set  up 
the  airfoil,  wake  and  farf  ioi  X5  turn'  takin9  account  of 

steady  solution  ?h?f  resets  ?n^Und;ry/0ndition8*  In  the 

for  the  column  of^V^h^e" 
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Gaussian ^elimination?  Ue^In^  ^ 

Uon^rr^ 

equation  is  linear  so  that  no  coluL  Ue^aUon^^red. 
After  each  column  is  solved,  it  is  roinvo^  „„• 

of^e^^i**!^^  th°  l0^  "ature 

hyperbolic  points  \h.  i eIllptlc  Points  and  w % .75  for 

B Sr  EfesigS" 

culation  is  a function  of  z along  the  airfoil^  In^he11  Clr~ 
txUonsnfor  theSs?eaSyCand  unstefd^  ^ date'  "Y"***^  solu- 

sldesrofUthf  g?id1X  Sn°t hiChl4  ^arrcond^Uo^o^ffvf i0n 
root  a sy^etry  houndary  c£S?t fVj"* 

as'iriirS 

the  solution  process  15  Updatod  Periodically  during 

in  fhIh?«SOlUti°n  pfocess  summarized  above  has  worked  well 
in  the  few  cases  calculated  to  date.  Convergence  for  in 

described  previously  for  STRANS  and  UTRANS  in  Referent  2. 
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Figure  4,  Coordinate  Definition 


The  development  of  three-dimensional  subsonic  farfield 
approximations  for  the  steady  and  unsteady  perturbation  poten- 
j ti^ls  proceeds  in  the  same  manner  as  described  in  previous 

work  for  the  two-dimensional  flow.  As  before  the  method 
involves  the  approximation  of  various  integrals  over  the  wing 
and  wake  which  result  from  the  application  of  Green's  theorem 
to  the  appropriate  partial  differential  equation.  Klunker7 
has  used  the  method  to  develop  asymptotic  solutions  for  the 
three-dimensional  steady  flow  and  his  result  in  the  following 
form  is  used: 

! I 

1 


^ ff  (*  > Y f z) 


ZSPAN 


/ / 


* (£  * H) d£dn 


-ZSPAN 


4tt  (y 7 f ? 7 ) 


ZSPAN 

^7  11  + / y(r|)dn 

-ZSPAN 


ZSPAN 


y2+z2  -»• 


X ->  —00 


/ dn 

-ZSPAN  (z-n) 2+y2 


for  x + 


where  R = ( (x2+k (v2+7 2 m **  4./,-  * . 
tribution  and  T(n)  is 

steady  perturba tion^po ten tiaiafol lows ithSOlUtK°n  for  the  un" 
and  is  now  described  in  some  detail  to!  Tet^d  o£  Klunker 
the  unsteady  perturbation  potential' ’ (^ttoSw^iStSS 


LI<,1  K*xx  + 0^  + 0^  - 2i(J*  + kR* 


= (<t>°  <p  ) 
x vx'x 


and  the  use  of  wing  and  wakfbo  h t0  the  linear  operator  L 

shock  conditions  results  iJ  the  fd? ?onditions  and  weak 
for  0:  results  in  the  following  integral  equation 
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(16) 


*% 


4>  (x # y # z ) - 


A<|>  (S  # n)  <1^  d£dn 


// 

wing 

/ Y<,,)/ 


+ 

SPAN 


ip  -iktf-l) 

S e dCdr, 


wake  integral 


Jff  ^ C d^dr|d^ 


where  ^ is  the  source  solution  to  L[4>]  = 0: 


• r 

'Mx,y,z;£,C,n)  = e 1 [ k (x‘^)  " 7=F  R 


4ttR 


JL 

/K- 


(17) 


where 


Vs1  (r  + k) 


R = 


S)2+k  [(y-rj  2+  (z-n) 2] 


The  use  of  the  source  function  <p  in  Equation  16,  neglecting  the 
volume  integral  as  a higher  order  term,  and  after  considerable 
manipulation  and  approximation  (as  x2+y2+z2-*®)  of  the  various 
integrals  results  in  the  following  farfield  solution: 
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in  the  equation 

theaendnof  can  be  evaluat  „ 

jt  ^ SGCti™-  lntegrand  and  is  p^nled'eT 

£*?  “«»  -stead,  farf  ■ , 

EE&1*  cs^the°nUmericai 

or  d 
The  use  of  b°-darles.  he  steady 

the 

i f results  in 


-Kf 


U+U? 


-ikr/uj 


+ ikr 


v=0  Cv+ikr 


-(cv+i kr)/u,| 


‘ (-  iSr)- 


v/ThT2 


-ikru  |lu,| 


«...  ..  A ® 

v=0  Cv‘,ikr 


where 


v -(cv+ikr)u  | /u  I 


f U+Y)6M2J  >/J 


t 

u = -1 
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4.0  PROG.'vAjh  descriptions 
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Computer  programs  TDSTRN  and  TDUTRN,  used  in  conjunc- 
tion, implement  the  theory  and  numerical  solution  procedure 
for  unsteady  transonic  flow  described  in  the  previous  two 
sections.  As  described  above,  the  boundary  value  problems 
for  the  steady  perturbation  potential  (Equation  9)  and  the 
unsteady  perturbation  potential  (Equation  10)  are  solved  in 
TDSTRN  and  TDUTRN  respectively  using  a finite  difference 
relaxation  procedure.  The  use  and  manipulation  of  magnetic 
tapes  forms  an  integral  part  of  the  operation  of  each  pro- 
gram as  well  as  serving  as  the  necessary  "data  link"  between 
the  two  programs.  As  a result  the  user  is  assumed  to  have 
some  familiarity  with  the  use  of  tapes  and  their  manipulation 
with  control  cards.  The  reading  and  writing  of  data  files  on 
magnetic  tape  is  described  in  the  next  section  and  motivated 
in  Section  6.0.  In  this  section,  the  logical  flow  of  the 
TDSTRN  and  TDUTRN  programs  is  described  and  a brief  summary 
of  each  subroutine  is  presented.  Both  programs  are  quite 
similar  in  logical  approach  and  operation,  so  that  they  are 
described  together.  Differences  between  the  programs  are  high- 
lighted with  appropriate  comments  as  needed. 


The  logical  flow  of  the  TDSTRN  and  TDUTRN  programs  are 
almost  identical  with  minor  exceptions  noted  in  the  description 
below.  The  calculation  is  begun  by  reading  card  input  and,  if 
a restart  is  being  performed,  a tape  dump.  In  TDUTRN  the  tape 
dump  of  the  steady  solution  being  perturbed  is  also  read.  . 
finite  difference  coefficients  and  airfoil  boundary  conditions 
are  initialized  in  a call  to  INITAL  and  subsonic  farfield 
quantities  are  initialized  in  a call  to  FARFLD.  If  a restar 
is  not  being  performed,  initial  values  for  at  all  grid  points 
are  determined  by  the  linearized  subsonic  or  supersonic  solu- 
tion. The  computational  cycle  is  executed  by  setting  up  the 
tridiagonal  equations  for  a column  of  grid  points  using  e 
mixed  differencing  finite  difference  equations.  The EJ 
are  solved  iteratively  in  TDSTRN  and  in  one  pass  in  TDUTRN,  by 
Gaussian  elimination  in  a call  to  TRI.  Each  column  is  solved 
and  relaxed  in  turn  proceding  through  the  grid  from  lef*  , 
right.  The  grid  is  swept  iteratively  in  this  manner  until  the 
change  in  $ for  all  grid  points  is  less  than  EPSGRD(l).  A call 
to  PRINT  prints  out  the  airfoil  pressure  coefficients  every 
NPRINT  iterations,  and  the  farfield  is  updated  every  NGFF 
iterations,  in  FARFLD. 


I 
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When  the  converged  solution  is  obtained,  a tape  dump 
of  all  relevant  input  and  calculated  quantities  is  performed 
and  a call  to  FPRINT  calculates  and  prints  out  the  airfoil 
pressure  and  force  coefficients.  Various  diagnostic  prints 
are  also  performed  in  TDSTRN  and  TDUTRN  after  every  grid 
iteration,  when  the  farfield  is  updated,  when  the  grid  is 
refined  and  when  a tape  dump  is  performed.  The  iterative 
procedure  may  also  be  terminated  when  the  maximum  number  of 
iterations  (NGRID)  has  been  exceeded.  In  either  case,  a 
final  tape  dump  and  final  print  are  executed. 

A summary  of  each  subroutine  is  now  presented. 


TDSTRN/TDUTRN 

These  are  the  driver  routines  for  the  respective  pro- 
grams. The  logical  flow  of  the  mixed  differencing  relaxation 
procedure  as  just  described  is  controlled  by  these  routines 
and  all  operations  including  input,  initialization,  finite 
difference  solution  and  output  are  performed  either  inter- 
nally or  by  calls  to  the  various  subroutines  described  below. 


D0UBLE 

(Not  in  present  version) . 


FARFLD 

The  subsonic  farfield  is  calculated  and  updated  in  this 
routine  using  the  asymptotic  solutions  for  the  steady  or  un- 
steady perturbation  potentials. 


FLP  (in  TDSTRN  only) 

This  is  a function  statement  which  contains  the  airfoil 
lower  surface  slope  distribution  used  in  the  linearized  tan- 
gency  boundary  condition.  This  function  is  called  from  sub- 
routine INITAL  and  its  value  at  each  grid  point  on  the  lower 
surface  of  the  airfoil  is  stored  in  the  FPL  array. 
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FPRINT 


This  routine  produces  the  final  print  and  is  called 
when  the  solution  has  converged  to  the  desired  accuracy  or 
when  the  problem  is  terminated  for  reaching  the  maximum  num- 
ber of  grid  iterations  allowed  (NGRID) . The  unsealed  pres- 
sure coefficients  above  and  below  the  airfoil  at  various 
specified  spanwise  stations  and  the  airfoil  force  coefficients 
are  also  calculated  and  printed  out  in  this  routine. 


FPU  (in  TDSTRN  only) 


This  is  a function  statement  which  contains  the  airfoil 
upper  surface  slope  distribution  used  in  the  linearized  tan- 
gency  boundary  condition.  This  function  is  called  from  sub- 
routine INITAL  and  its  value  at  each  grid  point  on  the  upper 
surface  of  the  airfoil  is  stored  in  the  FPU  array.  The 
doublet  strength  due  to  airfoil  thickness  (D0UB)  must  also 
be  given  in  this  subroutine.  This  quantity  is  defined  by  an 
integra)  of  the  airfoil  thickness  distribution  function 
(normalized  to  airfoil  thickness) : 


D0UB 


+ZSPAN  1 

■/  / 

-ZSPAN  0 


t (£, n)d^dn 


GAMFUN 


This  routine  performs  the  relaxation  to  update  farfield 
circulation  (GAMFF) . 


INITAL 


The  finite  difference  coefficients  AXl,  AX2 , BX1,  BX2, 

CX,  AY1,  AY2,  AZ1,  AZ2,  AX(DX),  AY (DY)  and  AZ(DZ)  are  computed 
in  this  subroutine.  The  airfoil  boundary  conditions  FPU  and 
FPL  are  also  set  here,  using  functions  FUP  and  FLP  respectively. 


PRINT 


This  routine  computes  and  prints  the  scaled  pressure 
coefficients  above  and  below  the  airfoil  every  NPRINT  grid 
iterations. 


TRI 


This  routine  solves  a system  of  tridiagonal  equations 
using  Gaussian  elimination. 


W^KE 

This  routine  solves  an  integral  used  in  the  unsteady 
farfield  solution  based  on  a rational  approximation  for  the 
integrand. 


1 


5.0  INPUT  AND  OUTPUT 


TDUTRN , AandSthePresul tinghQutput e^UihGd  t0  rUn  TDSTRN  and 
m this  section.  All  card  °f  ?h  pro9ran>  is  presented 

CDC  NAMELIST  package  with  a U^L^rd^"^ 


5.1 


TDSTRN  Input 


Recommended  and/or°typioalNvaluesWfo0nSidere.  in  three  sets, 
iables  which  control  the  nu^Hri,  wSOme  of  th,!  inPut  v«- 

ssss-1!.-: — 

•P—Sifc%rss5t0  rr^i^:fcSn? input 

n _•  


First  Set 


1 through  80  ' (Forma t^AlO^^’HThis3^  i[jformation  in  columns 
case  being  run  and  is  orin^  ! Can  ?e  USed  to  ^fine  the 
which  presents  the  final  conve?ged°resSl ts^  Pi>9e  °*  °Utput 


Second  Set 


rEois  Leame 

the  use  of  this  option  are  giveTTX^nSVKS'SXiS! . 


NAME 


ITAPE 


is^being^restartel  fl^"2  "--“^problem 
to  be  rlad  Irofa  d™pmtfpP:eV1°US  rU"  Whi°h  is 
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Third  Set 


The  third  set  of  data  is  read  in  under  NAMELIST  name 
$IN,  and  includes  all  of  the  variables  required  to  define  a 
problem,  and  control  the  numerical  iteration  procedure. 


NAME 

X 

Y 

Z 

IM 

JM 

KM 

ILE 

ITE 

JW 

KSPAN 

ZSPAN 

M8 

GAM 

DEL 

ALPHA 

GAMFF 

NGFF 


DESCRIPTION  

An  Array  containing  the  streamwise  grid  coordi- 
nates; IM  of  them 

An  array  containing  the  normal  grid  coordinates; 
JM  of  them 

An  array  containing  the  spanwise  grid  coordi- 
nates; KM  of  them 

Number  of  grid  points  in  the  streamwise  direc- 
tion (maximum  of  40) 

Number  of  grid  points  in  the  normal  direction 
(maximum  of  40) 

Number  of  grid  points  in  the  spanwise  direction; 
(maximum  of  20) 

I location  of  airfoil  leading  edge  (X(ILE)) 

I location  of  airfoil  trailing  edge  (X(ITE)) 

J location  of  airfoil  (Y(JW)) 

K location  of  wing  tip 

Wing  semi-span;  Z location  of  wing  tip 

Freestream  Mach  number 

Y,  ratio  of  specific  heats 

Airfoil  thickness  ratio  in  percent 

Airfoil  angle  of  attack  in  radians 

Initial  guess  for  the  spanwise  distribution  of 
airfoil  circulation;  to  be  used  in  the  initiali- 
zation of  the  farfield;  KSPAN  values. 

Every  NGFF  grid  iterations  the  farfield  is  up- 
dated ('''10) . 
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NAME 


DESCRIPTION 


0MEOAH 


0MEGAE 


0MEGAP 


EPSC0L 


NC0L 


EPSGRD 


KEPS 


NGRID 


NDUMP 


NPRINT 


Relaxation  parameter  for  hyperbolic  grid  points 
(^.75) 

Relaxation  parameter  for  elliptic  grid  points 
(^1.7) 

Relaxation  parameter  for  parabolic  grid  points 
(^.75) 

Convergence  criteria  for  column  solution.  The 
change  in  4>°  during  a column  iteration  at 
every  point  in  the  column  must  be  less  than 
EPSC0L  for  convergence  to  occur  (^5  x 10  5) 

Maximum  number  of  column  iterations  allowed. 

Note  that  if  NC0L  iterations  is  reached  without 
convergence,  a printout  of  the  degree  of  con- 
vergence is  given  and  the  calculations  proceed 
as  if  converqence  had  occurred  (^10) 

An  array  containing  criteria  to  control  grid 
convergence.  The  change  in  $°  at  every  grid 
point  during  one  grid  sweep  must  be  less  than 
EPSGRD (1)  for  convergence  to  occur. 

Set  equal  to  1.  (Not  used  in  current  version) 

Maximum  number  of  grid  iterations  allowed.  When 
the  number  of  grid  iterations  equals  NGRID  the 
calculation  is  terminated  and  a final  print 
given. 

Binary  tape  dump  frequency.  Every  NDUMP  grid 
iterations  current  values  of  all  variables  will 
be  dumped  on  tape.  Note  that  a tape  dump  occurs 
automatically  whenever  the  grid  converges  or  the 
number  of  grid  iterations  equals  NGRID  (set 
equal  to  large  number  if  a dump  of  only  the 
final  iteration  is  desired) . 

Every  NPRINT  grid  iterations  the  scaled  pressure 
coefficient  above  and  below  the  airfoil  is 
printed. 
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NAME 


DESCRIPTION 


IK  Setting  IK  * 1,  flags  the  use  of  a previous 

solution  as  an  initial  guess  for  the  current 
problem  where  the  Mach  number,  airfoil  thick- 
ness or  shape  and/or  angle  of  attack  may  be 
different. 

NKPRT  Number  of  spanwise  sections  for  which  pressure 

coefficient  data  is  printed  in  final  print. 

KPRT  K location  of  spanwise  sections  for  which  pres- 

sure coefficient  data  is  printed  in  final 
print  (maximum  of  20)  . 

ZE  Spanwise  locations  for  numerical  integration 

along  span  used  in  farfield  calculations 
(maximum  of  25).  This  permits  the  specifica- 
tion of  more  spanwise  points  than  available 
in  grid  (KSPAN ) to  increase  accuracy  of  the 
numerical  evluation  of  wing  integrals. 

NZE  Number  of  spanwise  locations  for  wing  integra- 

tion. 


The  input  data  listed  above  are  necessary  to  initiate 
a calculation  for  which  no  previous  calculation  is  available. 
Most  calculations,  however,  are  performed  as  restarts  using 
data  which  has  been  stored  as  binary  files  on  the  restart 
tape  according  to  the  format  described  in  Section  5.2. 

This  use  of  the  restart  capability  is  an  inherent  aspect 
of  the  recommended  computational  procedure.  Some  brief 
comments  describing  the  initiation  of  a calculation  using  the 
restart  capability  are  pertinent  at  this  juncture. 

It  is  noted  that  the  restart  or  dump  tape  (TAPE7)  may 
be  manipulated  in  any  way  desired  using  the  appropriate  con- 
trol cards.  In  general  the  tape  will  contain  data  from  many 
runs,  stored  as  individual  binary  files.  For  restarting  the 
TDSTRN  program,  the  desired  file  from  the  restart  tape  (TAPE7) 
is  copied  to  a disc  file  (TAPE8) . The  user  is  reminded  to 
rewind  TAPES.  TAPE7  is  then  positioned  at  the  end  of  the  last 
file  on  the  tape  so  that  new  dumps  can  be  written  by  the  program 
without  losing  any  of  the  old  data.  The  first  two  sets  of  data 
are  than  input  with  ITAPE=1.  In  the  third  set  of  data  the 
following  control  variables  are  needed  as  input: 

0MEGAH,  0MEGAE,  0MEGAP , EPSC0L,  EPSGRD,  NDUMP , 

NC0L,  NGRID,  NGFF,  PGFF , KEPS,  NPRINT , NKPRT, 

KPRT,  ZE,  and  NZE. 
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The  remaining  input  variables  are  stored  on  the  restart  tape 
and  need  not  be  input  unless  the  restart  option  is  being  used 
to  run  a new  case.  If  a new  case  is  being  run,  IK  must  be 
set-t?,1  whJ-ch  allows  the  Mach  number,  airfoil  thickness  or 
airfoil  angie  of  attack  (M8,  DEL,  ALPHA)  to  be  changed.  If 
the  airfoil  angle  of  attack  and/or  the  flap  angle  are  changed 

a,nc*J"ess  for  the  farfield  circulation  (GAMFF)  can  and 
should  be  made. 


5.2  TDSTRN  Output 


The  output  from  TDSTRN  consists  of  three  parts:  (i)  a 

continuous  commentary  which  describes  the  progress  of  the 
iterative  solution  procedure,  (ii)  a final  print  summarizing 
interest  from  the  final  converged  solution,  and 
(in)  a binary  tape  dump  of  all  pertinent  input  and  calculated 
parameters. 

The  continuous  commentary  consists  of  various  print 
statements  executed  in  the  main  program  TDSTRN  or  the  sub- 
routine PRINT  which  describe  the  current  state  of  the  solution 
as  well  as  the  occurrence  of  various  "milestones"  in  the  itera- 
tion  process.  The  only  print  that  occurs  every  iteration  is 
the  value  of  the  maximum  change  in  <t>°  throughout  the  grid 
during  one  grid  iteration.  When  a column  iteration  fails  to 
converge,  a print  occurs  which  defines  the  degree  of  column 
convergence  and  the  j and  k locations  of  the  most  poorly  con- 
verged point.  Every  NGFF  iterations,  the  subsonic  farfield 
is  updated  and  the  new  values  of  farfield  circulation  (GAMFF) 
and  airfoil  circulation  (GAMTE)  are  printed.  The  user  can 
examine  the  effect  of  degree  of  convergence  on  the  solution 
by  specifying  a print  of  the  scaled  pressure  coefficients  on 
the  upper  and  lower  airfoil  surfaces  every  NPRINT  iterations. 
Finally,  a descriptive  print  occurs  at  certain  milestone 
points  such  as  the  occurrence  of  a binary  tape  dump  and  solu- 
tion convergence. 

The  final  print  is  executed  in  subroutine  FPRINT  when 
the  solution  has  converged  to  the  desired  accuracy  or  when  the 
number  of  grid  iterations  equals  NGRID.  The  print  is  self- 
explanatory  and  includes  the  input  parameters  which  define 
the  problem  and  various  calculated  quantities  of  interest. 

The  calculated  quantities  are  of  course  based  on  the  final 
converged  solution.  The  section  lift  coefficients  are 
printed  out  as  well  as  the  upper  and  lower  surface  pressure 
coefficients  for  various  spanwise  coordinates. 
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The  most  important  form  of  TDSTRN  output  is  the  binary 
tape  dump  of  all  input  parameters  defining  the  problem  and  of 
the  most  recent  values  of  <f>°  at  all  grid  points.  A tape  dump 
occurs  automatically  if  the  solution  has  converged  to  the 
desired  accuracy  or  if  the  number  of  grid  iterations  equals 
NGRID.  The  user  may  also  specify  that  such  a dump  occur 
every  NDUMP  grid  iterations.  The  tape  so  generated,  not 
only  forms  a permanent  record  of  the  results  of  a calculation 
for  possible  future  editing  and  examination  but  also  forms  a 
necessary  part  of  the  computational  procedure.  Most  important 
is  its  use  as  required  input  for  a TDUTRN  calculation.  How- 
ever, it  may  also  be  used  to  restart  the  calculation  to  refine 
accuracy  or  convergence  or  be  used  as  the  initial  guess  for 

T throughout  the  grid  for  a similar  calculation,  as  described 
in  Section  5.1. 

The  format  used  for  writing  and  reading  the  binary  tape 
is  given  in  the  following  FORTRAN  statements: 


WRITE 

(7) 

NITERG , IM,  IM1 , JM,  JMl , KM,  KMl , JW, 

JWP 1 , JWM1 , ITE , ILE , KSPAN , KCAP , DEL , 
ALPHA , NDB , M8 , GAM, DYBU1 , DYBU? . DYBLl , 
DYBL2, D0UB, ZSPAN 

WRITE 

(7) 

(X(I)  ,DX(I)  , AX1  (I)  , AX2  (I ) , BXl  ( I ) , 
BX2 (I) , CX (I) , 1=1 , IM) 

WRITE 

(7) 

(Y  (I)  ,DY(I)  , AY1  (I)  , AY2  (I ) , 1=1 , JM) 

WRITE 

(7) 

(Z  (I)  , DZ  (I ) , AZ1  (I)  , AZ2  (I ) , 1 = 1 , KM) 

L=ITE*KM 

WRITE 

(7) 

(FPU (I) , FPL(I) ,PHIUB(I) ,I=1,L) 

WRITE 

(7) 

(GAMTE(I) , GAMFF ( I ) , 1=1 , KSPAN ) 

L=IM* JM*KM 

WRITE 

(7) 

(PHI (I) , 1=1, L) 

END  FILE  7 

Any  information  may  be  retrieved  from  the  tape  by  using  the 
appropriate  READ  statements  as  is  done  in  the  restart  option 
described  above. 
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TDUTRN  Input 


, fThe  inPut  for  TDUTRN  consists  of  normal  card  input  plus 
lnput  from  a binary  file  which  contains  the  steady  solution 

wh-bK  an  TDSTRN  run-  TDUTRN  also  has  a restart  capa- 

DrevinnQi'^H  1S  1Jp^e,!'ented  exactly  the  same  manner  as 

Spon  at  the  PndrnfSJh-n  SecJlon  5‘1  for  TDUTRN  and  elaborated 
upon  at  the  end  of  this  section.  The  required  input  is  now 

described  and  some  comments  are  presented  at  the  end  of  this 

AsCbeforeer^henin9^t<?  ^ read  °f  the  steady  solution. 

As  before,  the  card  input  is  described  in  three  sets. 

First  Set 

1 throughD8Si(Por^^8A!Sfai''in9  a"y  in£ormation  in 

Second  Set 

SCpNTRL^  SeCOnd  set  of  data  is  read  in  under  NAMELIST  name 


NAME 


n ape 


DESCRIPTION 

This  is  a flag  for  using  a restart  tape. 

ITAPE=0  means  the  problem  is  being  started  from 
scratch  (iteration  0),  ITAPE=1  means  the  pro- 
blem is  being  restarted  from  a previous  run 
using  the  restart  tape.  Note  that  a tape  is 
also  used  for  the  input  of  steady  results  inde- 
pendent of  the  value  of  ITAPE. 


Third  Set 


The  third  set  of  data  is  read  in  under  NAMELIST  name 


NAME 


DESCRIPTION 

An  array  containing  the  streamwise  grid  coordi- 
nates; IM  of  them. 

An  array  containing  the  normal  grid  coordinates* 
JM  of  them. 


An  array  containing  the  spanwise  grid  coordinate 
KM  of  them 
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name 


DESCRIPTION 


IM 

JM 

KM 

ILE 

ITE 

JW 

SMALLK 

KSPAN 

GAMFF 

NGFF 

PGFF 

0MEGAH 

0MEOAE 

0MEGAP 

EPSGRD 


KEPS 


uiiiw  4 


or  grid  points 
clon  (maximum  of  40). 

(maximu^of” 0)  f°lntS  in  the  normal  direction 
(maximum  of *20)  ^°ints  in  the  spanwise  directioi 

1 l0Cati°n  °f  airfoil  lading  edge  (X(ILE)). 

tion  of  airfoil  trailing  edge  (X(ITE)). 

J lo=ation  of  airfoil  (y (jW) , . 

Reduced  frequency  based  on  chord  = uc/u. 

K location  of  wing  tip 

ir^the'intial  i St  ioTo?  ‘ th'e ' f "“elation  used 

that  GAMFF  is  . cXiex  ^bef"fleld-  Note 

iation^n^the'farfield'is8  *5*  airf°H  '«cu- 

causes  the  farfield'to  b I SfiSS- (*Iof  alS° 

Relaxation  parameter  for  hyperbolic  grid  points 
Relaxation  parameter  for  elliptic  grid  points 
Relaxation  parameter  for  parabolic  grid  points 

convergencl?tamenIhMgeeina  co"tr°l  grid 

point  durinq  one  arin  ln  ^ at  every  grid 

EPSGRD(1»  convergenc^to^ccur?  leSS  tha" 

qual  to  1.  (Mot  used  in  current  version) 
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NAME 


DESCRIPTION 


NGRID 


NDUMP 


NPRINT 


IK 


XP 


Maximum  number  of  grid  iterations  allowed. 

When  the  number  of  grid  iterations  equals 
NGRID  the  calculation  is  terminated. 

Binary  tape  dump  frequency.  Every  NDUMP  grid 
iterations  current  values  of  all  variables 
will  be  dumped  on  tape.  Note  that  a tape 
dump  occurs  automatically  whenever  the  grid 
converges  or  the  number  of  grid  iterations 
equals  NGRID.  (Set  equal  to  large  number  if 
a dump  of  only  the  final  iteration  is  desired.) 

Every  NPRINT  grid  iterations  the  scaled  upper 
and  lower  surface  pressure  coefficient  per  unit 
angle  of  oscillation  is  printed. 

Setting  IK=1  allows  the  user  to  use  a previous 
solution  as  an  initial  guess  for  the  current 
problem  where  the  reduced  frequency  and/or 
mode  of  oscillation  is  different. 

Steamwise  location  of  pitch  point  for  pitching 
oscillation. 


ITYPE 


I0PT 


NKPRT 

KPRT 


ZE 


NZE 


Unsteady  mode  of  rigid  body  oscillation 

ITYPE=1  -*•  Pitch  about  XP 
ITYPE=3  -*•  Uniform  plunge 

Unsteady  formulation  option;  I0PT=O  for  low 
frequency  approximation,  I0PT=1  for  general 
frequency  theory. 

Number  of  spanwise  sections  for  which  pressure 
coefficient  data  is  printed  in  final  print. 

K location  of  spanwise  sections  for  which  pres- 
sure coefficient  data  is  printed  in  final  print 
(maximum  of  20) . 

Spanwise  locations  for  numerical  integration 
along  span  used  in  farfield  calculation  (maxi- 
mum of  20) . This  permits  the  specification  of 
more  spanwise  points  than  available  in  grid 
(KSPAN)  to  increase  accuracy  of  the  numerical 
integration  of  wing  integrals. 

Number  of  spanwise  locations  for  wing  integra- 
tion. 
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perturbation  potential  fc^at  the  solution  for  th*> 
solution  of  the  st^1'  imPlemented  in  TDUTRN  h un?teady 

is  accomplished  by generate"™;  K*3 

tape  generated  byTDSTRii  9- the  aPPropriate  fiie  Thls 

JP  the  restart  option^1"  m“ch  the  same  way  as  i ^Ump 

solutionf^  ineludingStheS tape ^ e a d *"o f ^ he t e a ^ ^ 

mesi“eereqSEV^DU™"  ^“se"  ££VU9hHy  "ore 

copied  from  TAPE7 eo"1"? • th®  desired  unstead^?  fUe'  TAPE8- 
are  rewound  fnr  fc°  a disc  file,  tappq  f tape  dumP  is 

tioned  at  thf  en^of  ^he^  thG  pro<*r™-  ' TAPE7^f8t!!nd  TAPE9 
for  acceptin a a ° ot  the  last  fiie  on  tho  C t7  then  posi- 

are  input  as  before  fKPG  dump‘  The  first  two6.,1!?  preParation 

set  o£  data  - sSSIS'Z 

5®«^SJ»*fPSSr  "S'  KGRID' 

^FT'  ZE,  and  NZE.  "KPRT,  KPRT,  I TYPE , 

Again  there  is  an  • 

and  I?^rd  f^\U!t"hCyl°"VorK1LW^dCeh  to  chanqe 

CUlatl°n  <GAMFF)  shou[dCbrmfdeeV'  9UEEE  for  theafa1r?ie(?“clr- 


TDUTRN  Qnt-p,,- 


tale^r3  • «S£2,5LSjSr  ££&&»*  or  TDSTRN 

print  and  binary 

™™'arr?0°pjyxd^f«en«  ?f  ^ 

Printed  out  in  fhl  ° ohat  the  real  and  varxables  in 

qe^eT^The^finai6—?”!^?  *the^unneeded^commen£ 

iur.^l0s.l;r'efn^  SSS  ssr- 

iasrsJtSss  sfs?  which  (S  - jf 

Also,  various  caaSula?LUnStead>’  solution  bei£o  “eady  s°lution 
calculated  quantities,  basedon^hfa^cin- 
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verged  solution,  are  printed.  These  include  the  real  and 
imaginary  parts  of  the  unsteady  contribution  (per  unit  angle 
of  oscillation)  to  the  aerodynamic  force  coefficients. 

Also  unsteady  contributions  to  the  upper  and  lower  surface 
pressure  coefficients  (per  unit  angle  of  oscillation)  are 
printed  for  every  computational  point  on  the  airfoil.  It 
is  again  noted  that  these  are  complex  so  that  the  real  and 
imaginary  parts  are  printed  out  in  order. 

The  other  form  of  TDUTRN  output  is  the  binary  tape 
dump  of  all  input  parameters  and  the  most  recent  values  of 
(Re^1,  Im^1)  at  all  grid  points.  As  before  the  tape  dump 
occurs  automatically  at  normal  program  termination  or  at 
the  users  discretion  every  NDUMP  iterations.  The  format 
used  for  writing  and  reading  the  binary  tape  is  given  in 
the  following  FORTRAN  statements; 


WRITE 

(7) 

NITERG  , IM,  IMl , JM,  JMl , KM , KMl , JWP I , 
JWM1,  ILE,  ITE,KSPAN,0MEG,SMALLK, 
DYBU1 , DYBU2  , DYBLl , DYBL2  , ND0UB  , XP 

WRITE 

(7) 

(X(I)  ,DX(I)  ,AX1(I)  , AX2  (I)  ,BX1  (I)  , 
BX2 (I) ,CX  (I) , 1=1 , IM) 

WRITE 

(7) 

(Y ( I ) , DY  ( I ) , AY1 (I) , AY2 (I) , 1=1 , JM) 

WRITE 

(7) 

(Z (I)  ,DZ (I) , AZl (I) , AZ2 (I) ,1=1, KM) 

L=ITE* 

KM 

WRITE 

(7) 

(FPU ( I ) , FPL  (I ) ,PHIUB(I) ,I=1,L) 

WRITE 

(7) 

(GAMTE  (I)  , GAMFF  ( I ) , 1=1 , KSPAN ) 

L=IM*JM*KM 

WRITE 

(7) 

(PHI (I) , 1=1, L) 

END  FILE  7 
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6.0 


PROGRAM  USAGE 


The  general  structure  and  usage  of  the  three-dimen- 
sional programs  TDSTRN  and  TDUTRN  are  very  similar  to  that 
the  original  two-dimensional  versions  STRANS  and  UTRANS. 

case'.i!:  fs  recommended  for  economy  sake  that 
user  initially  become  acquainted  with  those 
programs.  The  programs  are  documented  in  detail  in  Refer- 
ence 2 so  that  the  comments  concerning  program  usage  in  this 
manual  are  kept  necessarily  brief.  y 

al.  In  fcl}eir  present  configuration,  both  TDSTRN  and  TDUTRN 

niibPranfaX1I?'!Jm1°f  11  f 500  computational  grid  points  and  the 
number  of  grid  lines  in  the  streamwise,  normal  and  spanwise 
directions  must  each  be  less  than  40,  40,  20  respectively. 

g!?ratl°n'  TDSTRN  reciuires  70.5eK  words  to  load 
and  57.0 bK  words  to  execute  and  TDUTRN  requires  161. 7BK 

ti^nSw^°  1had  and  1L°*°«K  words  to  execute.  This  configura- 
tion was  chosen  so  that  each  program  could  fit  into  small  core 

ahlpa?eH°f  a/D<;  7600  comPuter-  If  greater  storage  is  avail- 
able and  used,  (Ex.  CDC  6600)  it  is  a relatively  simple  matter 

z°  FPU^FPL  tetcarray  Si26S  °f  the  primary  variables  PHI,  X,Y, 

Detailed  comments  and  suggestions  are  given  in  Refer- 
l c2ncefnin9.grld  design,  farfield  location  and  update, 
Tho  f °f  relaxatlon  factors  and  accuracy  and  convergence. 
nroSr*™ame  S01™*?*3  aPPly  to  present  three-dimensional 
n / d ^i11  n0t  bG  rePeated  here.  The  sample  cases 
presented  in  the  next  section  should  provide  some  guidance 
with  respect  to  such  items.  y 


7.0 


SAMPLE  CASES 


Detailed  input  and  sample  output  for  sequences  of 
TDSTRN  and  TDUTRN  runs  are  presented  in  this  section. 


7 . 1 TDSTRN  Test  Case 


A sequence  of  computer  runs  are  described  in  this 
section,  which  calculate  the  steady  transonic  flow  over  a 
6 percent  thick,  symmetric  circular  arc,  rectangular  plan- 
form  wing  with  aspect  ratio  8,  at  Mm  = .86,  a = 0.  The 
individual  runs  required  to  complete  the  calculation  are 
described  in  the  run  log  given  in  Table  1.  The  table  lists 
the  restart  tape  read  by  each  run,  total  grid  iterations, 
convergence  achieved  and  the  tape  dump  generated.  The  grid 
used  consisted  of  approximately  11000  points  with  IM=30 
over  -3.2  < x < 3.4,  JM*19  over  -5.4  < y < 5.4  and  KM=19 
over  0 < Z £ 6.0.  In  the  x direction,  16  grid  lines  were 
distributed~along  the  airfoil  chord  with  AX^  .06  and  in  the 
Z direction  10  grid  lines  were  distributed  over  the  span 
with  AZ'v*  .2.  The  runs  shown  in  the  log  implement  a "boot- 
strapping" technique  by  which  the  calculation  is  initiated 
at  a low  sub-critical  Mach  number  and  the  Mach  number  raised 
in  later  runs  to  the  final  desired  value.  _The  final  run  was 
taken  to  a convergence  of  A4>max  * 3.7  x 10  5 . All  runs  were 
completed  in  a total  time  of  65  seconds  on  a CDC  7600  which 
indicates  a computer  time  requirement  of  3.  x 10  CPU  sec/ 
grid  point/iteration.  The  final  convergence  achieved  is 
believed  to  be  more  than  sufficient  for  engineering  accuracy. 


Run 

M« 

00 

Restart 

Tape 

Used 

Grid 

Iterations 

Convergence 

Achieved 

Tape  Dump 
Generated 

IS 

n 

i 

i 

38 

- 3 

10 

IS 

2S 

.8 

IS 

19 

- 3 

10 

2S 

3S 

.86 

2S 

50 

j 1.9  x 10-" 

3S 

4S 

.86 

3S 

50 

3.7  x 10"5 

4S 

Table  1.  Sequence  of  Runs  for  TDSTRN  Sample  Case 
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7,1,1  lQPut  for  TDSTRN  Sample 

above  is  given  in^his Section.  °f  ^ TDSTRN  runs  described 
• Run  IS:  no  tape  read,  generate  file  is 

*** 3d  circular  arc*** 

SC0NTRL 
ITAPE=0 , 

$END 

$IN 

-9,i..i.i4,i.34.i.62;2;6lf^;!i;|f5ej'-76'-82. 

Z<1) 

IM=30, 

JM=19 , 

KM=19 , 

ILE=9, 

ITE=24 , 

JW=10, 

kspan=io, 

ZSPAN=2. , 

M8=.7, 

GAM=1 . 4 , 

DEL=.06, 

alpha=o. 0, 

GAMFF (1) =10*0. , 

0MEGAH=. 75, 

0MEGAE=1 . 7 , 

0MEGAP= .75, 

EPSC0L=5.E-5, 

EPSGRD ( 1 ) =1 . E- 3 , 

NDUMP=2000, 

NC0L=1O , 

NGRID=50, 

NGFF=2000, 

PGFF=1 . 5 , 

KEPS=1, 

IK=0, 

NPRINT=5, 

nkprt=io, 

xEa‘io:^?:3'4'5'6’7'8’9'10- 

NZE=2  , 

$END 
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• Run  2S: 


read  file  Is,  generate  file  2S 
***3D  CIRCULAR  ARC*** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH=. 75 , 
0MEGAE-1.7, 

0MEGAP=. 75, 

EPSC0L=5 . E-5 , 
EPSGRD=1 . E-3 , 
NDUMP=2000 , 

NC0L=1O , 

NGRID=50, 

NGFF=2000 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=5 , 

NKPRT=10, 

KPRT (1) =1,2, 3, 4, 5, 6 
ZE(1)=0. 0,2.0, 
NZE=2, 

IK-1, 

M8=0 . 8 , 

$END 


7,8,9,10, 


• Run  3S:  read  file  2S,  generate  file  3S 

**  * ** 3D  CIRCULAR  ARC*** 

$C0NTRL 

ITAPE=1, 

SEND 

$IN 

0MEGAH=. 75 , 

0MEGAE=1 . 7 , 

0MEGAP=. 75, 

EPSC0L=5 . E-5 , 

EPSGRD=l.E-4, 

NDUMP=2000 , 

NC0L=1O , 

NGRID=50, 

NGFF=2000 , 

PGFF=1 . 5 
KEPS=1 , 

NPRINT=5 , 

NKPRT=10, 
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KPRT(l) =1,2, 3,4,5,6,7,8,9,10, 
ZE(1)=0. 0,2.0, 

NZE=2 , 

IK-1, 

M8=. 86, 

$END 


• Run  4S : read  file  3S;  generate  file  4S 

***  3D  CIRCULAR  ARC*** 

$C0NTRL 
ITAPE=1 , 

SEND 

$IN 

0MEGAH= . 75 , 

0MEGAE=1 . 7 , 

0MEGAP= . 75 , 

EPSC0L=5 . E-5 , 

EPSGRD=l.E-3, 

NDUMP=2000, 

NC0L=1O , 

NGRID-50, 

NGFF=2000, 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=5, 

NKPRT-10, 

KPRT (1) =1,2, 3, 4, 5, 6, 7, 8, 9, 10, 

ZE=0. 0,2.0, 

NZE (1) =2 , 

IK=0 , 

SEND 


7.1.2  Sample  Output  for  TDSTRN  Test  Case 


The  following  pages  contain  a sample  of  the  continuous 
commentary  output  for  the  first  4 cycles  of  Run  IS  in  addition 
to  the  final  printed  page  of  all  runs.  Also  included  is  the 
complete  final  output  for  the  final  converged  result  (Run  4S) . 
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i.o  ci*cuc»*  ***C 


i ooon-oi  tirr  coi*mci»*t 


i 


"*•  **P  ClaCliL**  »ac 


7.2 


TDUTRN  Test  Cases 


A sequence  of  TDUTRN  runs  are  described  in  this  section 
which  calculate  the  unsteady  flow  perturbation  for  the  pre- 
viously described  circular  arc  rectangular  wing  oscillating 
in  pitch  about  the  leading  edge  at: 


M = . 86 


k = 0.0 
k = 0.1 


The  individual  runs  required  to  calculate  these  cases  are 
described  in  the  run  log  in  Table  2.  All  runs  used  the  steady 
solution  given  on  tape  dump  4S.  The  table  presents  the  reduced 
frequency,  restart  tape  read,  grid  iterations,  convergence 
achieved  and  tape  dump  generated.  The  runs  were  calculated 
using  the  same  grid  as  the  steady  runs.  They  were  performed 
in  the  order  shown  to  implement  the  "bootstrapping  technique 
for  getting  from  one  reduced  frequency  to  another.  The  input 
required  for  each  run  and  sample  output  are  presented  in  the 
following  section. 


Run 

k 

Restart 

Tape 

Used 

Grid 

Iteration 

Convergence 

Achieved 

Tape 

Dump 

Generated 

1U 

o 

• 

o 

50 

2.2  x 10-3 

1U 

2U 

o 

• 

o 

1U 

100 

1.9  x 10"4 

2U 

3U 

0.0 

2U 

50 

6.8  x 10-5 

3U 

4U* 

0.0 

3U 

86 

1.0  x 10'5 

4U 

5U 

— 

0.05 

4U 

50 

1.3  x 10~3 

5U 

6U* 

0.1 

5U 

200 

6.0  x 10"5 

6U 

Table  2.  Sequence  of  Runs  for  TDUTRN  Sample  Cases 
(•Denotes  converged  solution) 
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7.2.1  Input  for  TDUTRN  Test  Cases 


The  card  input  for  each  of  the  TDUTRN  runs  described 
above  is  given  in  this  section. 

• Run  1U:  read  file  4S,  no  restart  tape  read; 

generate  file  1U. 

***3D  CIRCULAR  ARC*** 

$C0NTRL 
ITAPE=0 , 

$END 

$IN 

X(l) =-3.2, -2. 2, -1.5 ,-1.02, -.67, -.42, -.24, -.1,0. ,.07, 
.14,  .21,  .28,  . 35, .42, .5, .55, .6, .65, .7, .76, .82, 

.9,1. ,1.14,1. 34,1.62,2.02,2.58,3.38, 

Y(l)  =-5.4, -3. 41, -2. 91,  -1.91,  -1.21, -.74,  -.43,  -.22,  -.08, 

0. , .08,  .22, .43, .74,1.21,1.91,2.91,3.41,4.3, 

Z(1)=0. ,.25, .5, .75,1.,  1.25, 1.5, 1.75, 1.9, 2.  ,2.1,2.25, 
2.45,2.75,3.2,3.85,4.75,6. ,6.8, 

IM=30 , 

JM=19, 

KM=19, 

ILE=9 , 

ITE=24 , 

JW=10, 

KSPAN*10, 

GAMFF (1)=10*(1.,0.), 

0MEGAH=  .75, 

0MEGAE-1 . 7 , 

0MEGAP=. 75, 

EPSGRD ( 1 ) =1 . E-4  , 

NDUMP=2000 , 

NGRID=50 , 

NGFF*=10, 

PGFF=1 . 5 , 

KEPS«1, 

NPRINT=5 , 

NNPRT*10, 

KPRT (1)*1,2,3,4,5,6,7,8,9,10, 

SMALLK=0 . 0 , 

IK=0 , 

XP=0.0, 

ITYPE=1 , 

I0PT=O, 

ZE(1)=0. , .125, .25, .375, .5,. 625, .75, .875,1. ,1.125,1.25, 

1.375,1.5,1.625,1.75,1.825,1.9,1.95,2., 

NZE=19, 

$END 
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• Run  2U:  read  file  4S,  restart  file  1U;  generate 

file  2U 

***3D  CIRCULAR  ARP*** 

$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH= . 75 , 

0MEGAE=1 . 7 , 

0MEGAP=.75, 

EPSGRD=l.E-4, 

NDUMP=2000 , 

NGRID=100, 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS=1 , 

NPRINT=5 , 

NKPRT=10 , 

KPRT= 1,2,3,4,5,6,7,8,9,10, 

ITYPE=1 , 

I0PT=O, 

ZE(1)=0. , .125, .25, .375  . 5 , . 625 , . 75 , . 875 , 1 . , 1. 125 , 1 . 25 , 
1.375,1.5,1.625,1.75,1.825,1.9,1.95,2. , 

NZE=19, 

IK=0 , 

$END 

• Run  3U:  read  file  4S,  restart  file  2U;  generate 

file  3U 

***  3D  CIRCULAR  ARC*** 

$C0NTRL 

ITAPE=1, 

$END 

$IN 

0MEGAH= . 75 , 

0MEGAE=1 . 7 , 

0MEGAP= . 75 , 

EPSGRD=1 . E-4 , 

NDUMP=2000, 

NGRID=50 , 

NGFF=10 , 

PGFF=1. 5, 

KEPS=1, 

NPRINT=5 , 

NKPRT=10, 

KPRT=1 ,2,3,4,5,6,7,8,9,10, 


ITYPE=1 , 

I0PT=O, 

ZE  (1)  =0 . , .25, .25, .375, .5, .625, .75, .875,1. ,1.125,1.25, 
1. 375 f 1.5,1.625,1.75,1.825,1.9,1.95,2. , 

NZE=19, 

IK=  0, 

$END 

Run  4U:  read  file  4S,  restart  file  3U;  generate 

file  4U 

***3D  CIRCULAR  ARC*** 

$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH= .75, 

0MEGAE=1. 7 , 

0MEGAP=.  75, 

EPSGRD=1 . E-5 , 

NDUMP=2000 , 

NGRID=100, 

NGFF=10, 

PGFF-1. 5, 

KEPS=1 , 

NPRINT=5 , 

NKPRT=10, 

KPRT=1,2,3,4,5,6,7,8,9,10, 

ITYPE*1 , 

I0PT=O , 

ZE  (1) =0 . , • 125, .25, . 375, .5, .625, .75, .875,1. ,1.125, 1.25, 
1.375,1.5,1.625,1.75,1.825,1.9,1.95,2. , 

NZE*19, 

IK=0, 

$END 

Run  5U:  read  file  4S,  restart  file  4U;  generate 

file  5U 

***3D  circular  arc*** 

$C0NTRL 

ITAPE=1, 

$END 

$IN 

0MEGAH=. 75 , 

0MEGAE=1. 7 , 

0MEGAP-. 75, 

EPSGRD=1 . E-5 , 

NDUMP=2000, 


NGRID=50 , 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS=1 / 

NPRINT=10 , 

NKPRT=10 , 

KP  RT=1 ,2,3,4,5,6,7,8,9,10, 

ITYPE=1, 

I0PT=O, 

ZE(1)=0. , .125,. 25, .375, .5, .625, .75, .875,1. ,1.125,1.25, 
1.375,1.5,1.625,1.75,1.825,1.9,1.95,2. , 

NZE=19 , 

IK-1, 

SMALLK=. 05, 

$END 

Run  6U:  read  file  4S,  restart  file  5U;  generate 

file  6U 

***3D  CIRCULAR  ARC*** 

$C0NTRL, 

ITAPE-1, 

$END 

$IN 

0MEGAH-.75, 

0MEGAE-1.7, 

0MEGAP=. 75, 

EPSGRD=1 . E-5 , 

NDUMP=2000, 

NGRID=200, 

NGFF=10, 

PGFF-1.5, 

KEPS—l , 

NPRINT-10, 

NKPRT=10, 

KPRT=1 , 2,3,4,5,6,7,8,9,10, 

ITYPE-1, 

I0PT-O , 

ZE(l)-0. , .125, . 25, .375, . 5 , . 625 , . 75, . 875 , 1 . , 1 . 125 , 1 . 25 , 
1.375,1.5,1.625,1.75,1.825,1.9,1.95,2. , 

NZE-19, 

IK-1, 

SMALLK- . 1 , 

$END 
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7,2,2  Sample  Output  for  TDUTRN  Test  Cases 


The  following  pages  contain  a sample  of  the  continuous 
°?tPU?  f°r  the  firSt  9 ^clesof  Run  10  in  adSi- 
the  fJnal  Printed  page  for  all  runs.  Also  included 
in  the  complete  final  output  for  the  final  run. 
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RUN  3U 
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RUN  4U 
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• FINAL  OUTPUT  FOR  RUN  6U 
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• FINAL  OUTPUT  FOR  RUN  6U  (CONT'D) 
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• FINAL  OUTPUT  FOR  RUN  6U  (CONT'D) 
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appendix  a 


r 


c 

c 


c 

c 


1 PRTAPf ft) TDSTRN  nNBuT»nuTPUT»T*pES»!NPHT,TAPt6»0UTPUT,TAPF7, 
PEAL  *CAP,m8,  Iwtk/G 

? KJo.rp.jtJooi.rPLjJoo?  phihn!«!a!  52* : 47?(?0>»*^«’).vuo)# 

1 ALPHA  t6AMf  F,p*f  GAH,OMrGAE,6“EGAP  f PAOHl  * "**  G<  * ' * 

,w* 

data  AL,8feT  /.5,,5/ 

FlLrE?IiE«#PSns??TAMC?rjF7MJTD?I?  [To  S”™!1,  [1°*  T*PE7  T°  4 DI9C 

ivsw ™ ?:Sr°:,?Su!?^KiE.s;  smt 

!*f!m  ‘TmE(!,,!.,|B) 

PEAD  (5, COnTRl ) 

«AiM;:Arfc5QrfjTJ2tT?,i? 

2 DOUB,  ZSPAN  L'*LPH4#NDB»M|,'Cah»O^U1,DVIIu?#DV0L1.DvBL?, 

EJ?e"’  C7cn.Dfc«,.Ai,S;)^S{fS,:Il:i;j;; 

read  (8)  crPii(i,,FPt«i),PMiuem.i.i.Lj 
T-lS*  j2lKi<;AMTE  f 15  'GANFF  C * 5 * J*1  *KSpAN5 

HEAD  (0)  (PH! tI)#I*1#L) 

I K*0 

READ  (5,IM> 

AlRrOlLnTHICKSIfSSESEAND/npnANCtSF!4nrTr5  DIFFr*ENT  h*CH  NUMBMs, 

THE  HEw  H8,  DEL  AhS/JJ^LPMA  GUS  ^ 9URF  YOU  HAVE  INPUT 

IE  (IK.EQ.O)  GO  TO  ! 

CALL  EARFlD 

1 CONTINUE 
8K*8QRT (KCAP) 

DO  2 J«1,KBPAN 
SAMTEl ( I )«GAMTE (T) 

2 CONTINUE 
WRITE  (6,900) 

WRITE  (0,90!)  NITERG 
NITERGpO 
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GO  TO  15 

C START  PROBLEM  from  SCRATCH 
10  CONTINUE 
READ  (5, IN) 

KC  APa ( j c«MS**2)/( ( 1 .♦GAM ) * DEL4MA**2)a* ,6666666667 
SKaSQRT(KCAR) 

00  5 i«i,kspan 

GAMTE(I)«GAMFF( J) 

gamtci  (n«r.AMFF(i) 

1 CONTINUE 
NITERG«0 

NDBaO 

IM1«IM-1 

kmi«km-i 

JwP  1 a J w + 1 
JWM1 BjW«l 

C INITIALIZE  FINITE  DIFFERENCE  COEFFICIENTS  AND  FARFICLD 
call  inital 

CALL  FARFLD 

c initial  guess  for  subsonic  case  (interior  only) 

00  20  Kai,KMi 
MP«IM*JM*(K-1 ) 

Z2«Z (K)a*2 
00  10  I«2,IM1 
M"MP4 ( I *1 ) a JM 
X20(I)*«? 

C0N»*X(I)*D0UB/(6.?MlB53) 

00  (10  Ja2,JNl 
L«M4j 

M9QRT  (X2+KCAP*  ( y ( J)**2+Z2) ) 

IF  (R.EQ.O.J  GO  TO  41 
RHI (L)*C0N/R**3 

IF  (ASS(PHI(L)),GT,1,)  PHI(L)«SIGN(|,,X(I) j 
GO  TO  40 
41  CONTINUE 

fhi(l)»phi(l-jhj 
40  CONTINUE 
30  CONTINUE 
20  CONTINUE 
L«ITE*KM 
00  5 I«1,L 
FHIUB ( I )>0 • 

5 CONTINUE 

H»( ILE»2)*JM4JW 
KK«(ILE-1)*KM 
00  47  K»1,KSPAN 
L«H4lM*JM4(K-i) 

FHIU8(KK4K)«PHI(|  ) 

47  CONTINUE 
15  CONTINUE 
WRITE  (4, IN) 

WRITE  (6,900) 

CPCPBaDEL** , 6666666667/ ( (1 , 4GAM) *MB**2) **,3333333333 
WRITE  (6,913)  KCAP,CPCFB 
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r 


RE. 

75 


FOR  GRID  ITERATION 


KGRD«1 
CYCLE  POINT 

continue 

£RROR*0  • 

rit«ni?erg 

NlTERG*NlrERG*l 

IE  (MOD(N!TERG,NPRlNT) ,EQ,0)  CALL  print(ntts 
IF  (HOD(Nm,G.,GFF)..C.O  GP  TO S” 

CALL  GAMFun 
CALL  EARFLO 

76  CONTINUE 

8ECImJh2?S«3h  THE  "L‘NE9  « ««CTI«.> 

DO  J00  K»i , km j 
MP*  I M JM* ( K *1 ) 

CHEXFOTL»0 AXRF0XL  CC0RWCCT  pLANE) 

mmJA  <k*LE.RRPAN)  ifqil.i 

BEGIN  LOOP  on  a GIVEN  PLANE  fx  DIRECTIONS 
00  200  I»2, x M 1 »IRtCTION) 

MBMP6(X«|)*jh 

CMECK  FOR  AIRFOIL  (CORRECT  COLUMN) 

IrLAG«0 

SAVE  THIS  COLUMN  OF  PHI 


201 


COLUMN  ITERATION 


396 

ITERATION 


on  201  J*2 # Jm 1 

L«M*J 

FHIOG( J)»PHI (L) 

CONTINUE 
NITERC«0 

LOOP  BACK  POINT  FOR 
250  CONTINUE 

NI*ERC«NITERC*! 

• AucF.i2I!tRC,5T,NCl5L)  Co  T0 
SAVE  PREVIOUS  PHJ  FOR  COLUMN 

00  202  J»2,JMl 

L«M4J 

FHIOC ( J)»PHI (L) 

202  CONTINUE 

0M  CnLUMN  DIRECTION) 
00  300  J»2,JM1 

CALCULATE  CELL  IN0ICE5 
L«M»J 
L*«L*JM 
LL»L»JM 
LLL»LL-JM 
IE  (I, EG, 2)  LLL*LL 
LA«LM 
LBpL-1 
LEpLMMJM 
L8K«L-IMJm 
XE  CK.EQ.l)  LBK.LF 
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i 


o o o 


I PHIPbPhI (LP) 

TPHIL«PHI (LL1 
TPHIU«PHI(LLL) 

TPMIBKbphi  fe.RK ) 

(Ifoil.fq.o.pp. j,ne,jn>  co  tp  in 

if  • • < c itf-i )*kh^k)«phz cll> J 

IF  (I.EQ.TTE*!)  PHI U LL)«.5*(PMIuBf ( ITC-?)*KM  + K)TPNT ft  i i 1 % 

i1FpMifiEii^Ii!f2iJn!:^!Tr,4ND*J,EQ*J“*4Nf)*K*EC*KSP*N4n 

1 P I(LBK)*,5*(PMHlH(n»l)*t<M4KSPAN)4PHI(LBKJ} 

: t :i  5o*J,7,,—>  • • 

TP  CO  TO  3M 

PAPT■0L4G,ER'1**Nn*J,E,^*Jk,M1,  Gn  Tn  370 
e M;D!?{2;,0',,rp,L>r’-'’’ co  T0  «* 

L p*9t“«5*(AV1(J)«/»Y2(J))*5IGI 

1»!  CONuSuf”  ’ 

IF  (VfJJ.LT.O,)  GO  TP  3R0 


^ FLLIPTIC  DIFFERENCING 

OMEGA  (.?)  BnMtGAF 

IF  (J.E0.2)  GO  TP  305 
IF  (J.EQ.JMi)  rp  TD  SOU 
GO  TO  390 
BOTTOM  BOUNDARY 
303  CONTINUE 

D(J)aO(vn-AY?(J)*PHi(in) 

GO  TO  390 
TOP  BOUNDARY 
30«  CONTINUE 

0<J>BD(J)*AYt (J)*PHI(IA) 

GO  TO  390 
320  CONTINUE 

IF  (Vm.GT.O.)  GO  TO  390 


- " > ■ > "”""PVVVVfft9tffffAARd 

C *****«*•«*•*  HYPERBOLIC  DIFFERENCING 


OMECAf jjbomEGAH 
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(J) 1 15  #fPH1  (L5"PMI  fLL))-CX  (!•?)*  (Phi  (LL).pHIfUL5) 
C(J)«/ky2(j]  f 1-1  5"*V1  f J5-4v2(J)-A?|  («)-A22fK) 

viiS 

l*  (J.ta.jMn  go  rn  123 

GO  TO  390 

C BOTTOM  BANDAS  V 
322  CONTINUE 

G(J)*0(J)«4yp(j)#PHI(, 0 
. CO  TO  390  u J 

c TOP  BOUNDAPV 
523  CONTINUE 

C **•«•••••••*«**. «*#,,#ll  ^4R4B0LIC  OIFftRERCINC  *•*•******•*•***•***••• 

390  CONTINUE  ************************************* 

OMEGA ( J)*om£GAP  *•**••*• 

A( J)«Avi ( j) 

IF  CJ.EQ.Jm^  Gn  TO  3413 
CO  TO  390  J 

" BOTTOM  BOUNDARY 
392  CONTINUE 

DCJ)«0CJT.Ay2(j).pH,aB. 

GO  TO  390  L 1 

C TOP  BOUNOAy 
391  CONTINUE 

0CJ)«0(  JJ.Avt  (J)*PMI((.*j 
GO  TO  390  L 1 

C ***********************»*****#(l# 

c ************************8H0^K  OIPPEPENCING  ****************** 

"*  5.~.n.ts.t  

AUJPAYl  (J) 

•AVI  (J)aAY?(J)« 

CCJ)«AV2(J) 

********** **,#,##  !******************************#j> 

soun^.r,  coNomS“*‘*“ ***** 

MO  continue  ****** 


i'  O'v.it.o”' cn^T^ ' ^ " j“p“ItL5>_*1<S(l)*CPMici) >phi cld > 
c eu{St!c<J>'1’t‘,‘1  00  70  555 

OHEGA  f J)»omEGAE 
A(J)«f5VBui 

551  CONTINUE 

: M 7°  5” 

0HEG4(J)«rjHFG*H 

c f 5-Oreui-AZj  CK).AZ?c«o 

PAPABOl ic 

552  CONTINUE 
OHEGA  f J)«n*«EGAP 
A ( J) «OvBUl 

c(j!ir#8i<,n"n*DvBu,*4z,<,(,*422(,<> 

SHOCK  POINT 
555  CONTINUE 

QHfGA  f J ) ■riMfcGAfc 
A(J)«0v»jui 

i D * ciS  ?*  * ^K{7JrJi:ss;{“» « 'Vi^y***  >-»n»  *tt , v f j 

l 572(K)*PNI(L*U))  5 f HI (LL )«PHI (LLL)))"(A?1(k)*phI(lF 

GO  TO  390 

».  csiiji, HHHSiE 

tLl!JTK,J,’lT‘0'’  '°  T0  573 
OHEG4(J)«OMtGAE 

•AZ2(K) 


on  TO  390 
*71  continue 

e mv»mboJ!«,'gt*®*’  50  Tn  171 

OMfG4( JJbomeGAh 

c PARABOLIC 

372  CONTINUE 

OMEGA (J)cQMEGAP 
A(J)*0, 

c SHOCK  POINT 
I7S  CONTINUE 

OMEGA  f J)*0ME6AE 
A ( J) *0 ( 

C(J)OV0U1’M'*VV#(BMCn40,'2(In4BtT#vfJJ*BM(I-l).A7l(K).4??(K) 

1 V ( J)#(B)f f ( i»  j j I ) *PM!  (LL) ) ABET* 

3*0  CONTINUE 
C BOO V BOUNOAPy  JmJk 
A ( J ) *0  , 

C(J)«0# 

0(J)«PHI(L) 

1B0  CONTINUE 

BMI(lh)«TPhi« 

*MI (LL)«TPM|l 
BHI(LLL)pTPHILL 
BMlaBKJ.TPMIBK 
300  CONTINUE 

: T",e0it‘cso?;i«;:;r ,s  5tT  “n-  *»*  »» **««, 

CHt55  m » W 

L»M*J 

JERROR«J 

EPRC«PHIOC(J).PHI(LJ 
jm  eo*Ti!!*e',"e,,ST,t',*e00  s0  t0  *50 

SB«  CONTINUE 

eMH.j8”S;K^,MS,tMos,::2  axwrjgss**" 
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IS  SET  NOW  SOLVE  FOR  COLUMN  OF  PHI 


i 


00  396  JS2.JM1 
L*M*  J 

CBB«0MEG*(  J)*(PMT  (L)-PHTfJG(J)) 

PHI  (L)«PHIOG( J)4fBB 

IP  (ABS(ERB),LT,AA3(ERRnR))  GO  TO  396 
EP90R«EHB 

LEBBOP«L 
196  CONTINUE 

IP  (IFLAG.NE.l)  GO  TO  200 
L«M*JW 

PHI  (L)«PMI  (L«lHC'Y(JwMn*(PMI(L-t)-PHl(L-2))/DY(Jw-?) 

phiub(n)«pmi  a^n-ov  ( jw)*(pwi(l^?)-phi(l*i)  )/oycjhpi) 

IP  (I.EO.TTE)  GAMTE(K)«PHIUB(N)-PHI(L) 

200  CONTINUE 

100  CONTINUE 

C PRINT  OUT  ERROR  AFTER  GRIO  SwtEP 

WRITE  (6,905)  NITERG, ERROR, LEBPOB 
IF  CARS(ERROR),LT.10,)  GO  TO  101 
WRITE  (6.915) 

STOP 

101  CONTINUE 
I0OUB«0 

IP  (AHS(EHROR).LE.FPSGRD(KGRP))  GO  TO  R0O 
IF  (NITERG.EQ.NGRIO)  go  to  RIO 

ip  cnoO(niterg.nduhp) ,cn,o)  go  to  rio 

GO  TO  75 
R00  CONTINUE 

KGRD«KGRD4l 
IOOUB«1 
GO.  TO  Rio 
R01  CONTINUE 

CALI  GAmPiin 

WRITE  (6.910)  NITERG.IWING.GAMTF (1).GAHFF(1),GAMTE(KSPAN), 

1  GAMFF (KJPAN) 

CALL  fppimt 
WRITE  (6,900) 

< WRITE  (6.906)  NITFRG 

call  onuBifc 

WRITE  ( 6 . 9 1 R ) IH.JM, Jw,KM,!LE,lTC,*<SPAN 
WRITE  (6.902) 

WRITE  (6.903)  (*(I),I«1,IM) 

WRIir  (6,911) 

WRITE  (6,903)  (V(I),I«1,JM) 

WRITE  (6.916) 

WRITE  (6,903)  (7(1), 1*1, KN) 

GO  TO  75 
•10  CONTINUE 
C TAPE  DUMP 

NPITE  (7)  NITERG,IN,IP1,JM,JM1,5M,KM1,JW,JWP1,JWM1,ITE,II!1 — - 

1 KSPAN.WCAP.OFL, ALPHA, nob, MB, gam. DYPU1,DVBU2,DYBL1,DYBL2. 

2 OOUB.ZSPAN 

write  (7)  (X(n,OX(n,Axim,AX2(n,BXl(!>,BX2(I),CX(n,I«l,!M) 
WRITE  (7)  (Y(n,0Y(n,AY|(!),AY2(n,!«l,JH) 

WRITE  (7)  (Z(I),BZ<n,AZim,AZ2Cn,!«l,KM> 

L*ITE*KM 
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il 


i 


r 


L«lM*JM*KM  e”’#GAMFF(  I),I«1,KSPAN) 

wTSu["/',lln'|,|,tl 

If**,1*  *6'9075  NJTFRG 
CAJ.I  PRINT(NITERG) 

I <«GBD.GT.KEPS)  go  TO  <120 
rF  fNITEPG.EO.NGPJO)  G0  ™ <*3a 
i?  (loouR.to.,,  G0  ™ “3' 

*28  CONTINUE 

JJITE  (6,908)  NITERG 
CO  TO  950 
*58  CONTINUE 

**ITE  (6,909) 

♦88  FORMAT  ( 1 H ] ) 

♦81  FORMAT 
♦82  FORMAT 

♦89  Fn^M4T  (10C*3»5J 

"i1",:', $\K;;v  “»•*"“*«”;  coium,<*,“*  '««•  «>  convc»ge. 

" ^o'oilA]  THE  H.EIHOH  E„0H  ..E„.S.  Qc 


^ITLRG 


> /»* 

#/,* 


REST“'Tt0  ” ite«*t:oh.„, 


UH  ,/ 


( 1 M 
(1M 


NUMMER  of  *ODE5  is  BEING  DOUBLED  AT  ITERATION* 

S-cSIBW^iB-uu, .?  lm 


i curbed 

♦86  FORMAT 

, * I5J 
*87  FORMAT 

♦88  FORMAT  m 

1ATI0N*I5) 

'MsAk  i5fr,SfS"jrSi""TS^*****T*0"* MM m*  ’rlCHCB^ 

*u “owU*;?;  

♦12  FORMAT  (8 A 1 0 ) * VfJ5»J«l»J"*) 

,,i»0*ciJeiJ»i(;/;IEfj‘,^*,,TV  p»«*»eteh 

l}2'22„. 

1 • ITE  •* 1 9 * KSPAN  .*19)  J 


CAS 


rAMrc  . * --"-I 

CAMFF ( | ) **Cl 3,5*  GAMTE(KSPAN 


(K) 


( 1 H , / . 


SOLUTION  is  DIVERGING 


'/f*  7(K),K«1,KM*) 


**tI3,5,/,*  SCALING  facto 
••19*  km  b*I9*  ILC 

0 THE  problem  is  bfing  termin 


♦ 15  FORMAT 
1*TED*) 

♦l*  FORMAT  (1H 
*58  CONTINUE 
CALL  FPrInt 
End 

subroutine  double 
”Al-  KCAP,M8,!WING 

j 

kb 

-71- 


RETURN 

END 

subroutine  FARFLD 
RIAL  KCAR, Mg, I m T n G 

1 s*°c«o?fcIf«o?rAvi  taoj  #0K* f » 

2 Z(  20)  •FPIJC800)  ,FPU  f 800  ) . FHIUB  f BOO  I 1 2 J lif  2 l*  X f a° J ' Y f a0  5 * 

common  /TNTe»p/%E;U  f;::eAMTCt*«>-”"-e*H'F(?0) 

; SUMONIC  hrpield 
' CALCULATE  rInG  INTEGRAL 
C0Nl«00UB/6,2fl3iF53 
IWING«0# 

00  10  I*2»  K3RAN 

to  MimJ?,NS*,,*t“"TC<,,*G*"Te(t*'»*6'ci-n 

JS:;^n:5t:;mGrTC(Ksp‘N,*(isp‘N-z(A^*">>> 

Z«Z(KM) 

MPljM* 

Z2«Z(km)**2 
00  20  I«1,IH1 

M«MP4>  f 1-1  )*J* 

Y2«X(I)**2 
CON3*«X(I)*CON1 
00  25  J«1,JM 
L«M*J 

V2«V(J)**p 

RaS0RT(X24KCAP*(vp473)) 

RHlTaC0Nl/R**3 

20  CONTINUE 
X«X<1) 

R2«XM  )**2 

CON3«-X(n*CON! 

00  30  K«1,KM1 
MR*IM*JM*(K-1 ) 

Z2«Z(K)**p 
00  35  J«1#JM 
L*RP*J 
Y2«V( J) **2 

R*SQRT(X2*KCAP*(y2*72)) 

FHIT«C0N3/R**3 

^ J 0 »» AND, Z2 , EG • 0 « ) GO  TO  54 

60IT0)35HIT4C0Na*V(J)*<1,t)<n5/ff,/(Y2*Z2J 

5*  CONTINUE 

Rhkljiphit 
M continue 

SO  CONTINUE 
Y»Y<1)  AND  Y*Y ( JH) 


i 


TO 


°o  ao  ro«i # ^ 

*[/!». £0.8)  on 
«■! 

50  TO  ap> 

•*  CONTINUE 
J»JM 

a*  CONTINUE 
°CI  as  K«1#«hi 

«■£? 

s !”:ess;ji/,/,mM» 

■ .ffiSHSSx:;"* 

•*  eSiftiS£M,T*eoN»*«* •♦««/») 

CONTINUE 
. a0  CONTINUE 
C X«X(IM) 

**»XflM)*** 

00  50  K«1,km 

JO  TO  I»1,WEE 
4fIi«ZECI) 

T®  CONTINUE 
NENOiNZE 
IELIP»0 

*2  CON^Jur5 
NEnO*K5PAn 

. IELIP.J 

T1  CONTINUE 

L»M4J 

JO«;vu, 

so  to  t; 

« ?0°,f}4<,,‘*‘EN0 

‘«»TK!|,,W,Z,,fE,>  60  TO  S3 
JO  TO  M 
*J  CONTINUE 
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n*(6AMTEClZ|). 


1 5*MTF(IZE-1 )))/((Z(K).A(KK))**2*v2) 

IF  (KK.EO.n  GO  TO  611 

AlNT«AlNT*,5*  (AnF**OLO)* ( A (KK )•« f KK>1 ) ) 

64  CONTINUE 
OlO«ANEw 

6 1 CONTINUE 

MNT«AlNT*,5*ANEW*(ZSPAN.Z(K3PANn 

65  CONTINUE 

»■ SORT  ( X2*KC AP*(V2* 72 ) ) 

PHlTs*v(lH)*C0Ni/H**3 
PmI(L)«PHIT*CON3*AINT 
60  CONTINUE 
50  CONTINUE 
RETURN 
END 

function  flp(xv,7z) 

PEAL  KCAP, MS, IWING 

»COir^i,»?fi:Ty-o,f(a0),Ov<afl,fO7(J!8,»**tt«05»FF2(«0),8*i(ao). 

1 9 <7  (6  0 ) , C I ( <l  0 ) , A V t ( 40  ) * A V2  ( 60  ) , A 7 1 ( 20  ) , AZ2  ( 20  ) VCaAi  vfartt 

J 7(PO),FPl»CP00>fFPLC«00).PMluB(800)jH,;!li;jJ®];f(J!  il!  52* 

] Jn»UJ-1.ITE.ILF,K3P4N.OvPu1,OvPu?;o;Jl1^ !dvBL2 X%mJ  Ul  i* 

« CAN. KCAP. NOP, TITLE  <*,  Do„0  ,w  nG,7SPAN.Umt  ipiloA?*  L*  # 
AIRFOIL  LO-ER  SURFACE  SLOPE  DISTRIBUTION * MT*' 

FLP»a.*txx-t5) 

RETURN 

END 

subroutine  FPRINT 
REAL  KCAP.M8, IWTNG 

COMMON  /OFl.T  A/  DV(40),Dy(40),OZ(20).AXl  (uni  iyiUiAi 

l cl!  «cir;op  XaSI  Vl:°v,uf#0m,'0VBL**4LM4*®p^^* 

COMMON  I Jit*  CJJ  •OOU»*  X“lRAC#/8FA4<#S|«^»T#«P»Tf 

/COEFF/  A(40),S(40),C(40), 0(40), PH 1(11500) 
r2r.2^.G*  6AMTt*<285iCAMTE (20),PGFF#gamFF(20) 

WRITE  (6,9001 

WRITE  (6,901)  (TITLF(I),I«1,S) 

WRITE  (6,902)  MS 
WRITE  (6,903)  KCAP 


WRITE  (6,901)  (TITLF(I),I«1,8) 

WRITE  (6,902)  MS 

WRITE  (6,903)  KCAP 

WRITE  (6,904)  DEL 

WRITE  (6,905)  ALPHA 

WRITE  (6,907)  7SPAN 

WRITE  (6,906)  CPCPB 

WRITE  (6,916)  CPCRIT 

WRITE  (6,912) 

WRITE  (6,915)  (X( I ) , I«ILE, IT€) 

cliftro, 

Cmomsq , 

00  10  W«1,KSPAN 
U««MP4ITE*JM*JW 

PHI (IJK)«PNI (IJK)*PART 
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(6,907)  7SPAN 
(6,906)  CPCPB 
(6,916)  CPCRIT 


(6,912) 

(6,915) 


(*(I),I«lLE,m> 


i 


LR*(ILE«2)*KM+k 
PHIUB(|.P)«phi  (m 
LP«ITE*KH4K 

W«MP*(I.1)*JM 

l»W+Jw 

LP*(I»1 )*km+k 

^cijHr^r  axi  cl^cpHTfL^JM>"pHI^L,^A¥2^I5,fc(pHIcL5•pwl<L-J'w,5^^ 

I-.5“r®T?*>  so  to  si 


1 


Co  TO  22 


y CK.GT.l) 

C(I)«A(I) 

OCIJ»0(I) 

21  Continue 

C2»C!*X(I) 

]p  fl.CT.ILE) 

CL«C1*X(ILE) 

CM«.S*C2*XflLE) 

CO  TO  23 

22  CONTINUE 

cMl£!;*,2*fCucl0)*£,x<J*i) 

« eST«JaS**e**e*0,*#*,f-t » 

cio«ci 
C20aC2 
*0  CONTINUE 

Jr Jj#?°!n  60  t°  u 

CH0MlcMOHri#?JfCL,fCL0,*0Z(KB*> 

CLO«CL 

CMQaCM 

JO  12  Na] # NKPRT 
J!mJ^RT(N)#NE»k)  CO  TO  12 

W»TT?TJf#IGAMTr(K,*CPCPB 
RITE  (6,906)  Z(K),GAMPRT 


WRITE  (*,913) 
WRITE  (6,915) 
WRITE  (6,914) 
WRITE  (6,915) 
00  TO  10 
12  CONTINUE 
10  CONTINUE 
WRITE  (6,9001 
WRITE  (6,901) 
WRITE  (6,902) 
WRITE  (6,903) 
WRITE  (6,904) 
WRITE  (6,905) 
WRITE  (6,907) 


(8(I),I»ILE,ITE) 

(A(I),IbIlE,ITE) 


<TlTLE(I),ial,8) 

KCAP 

OEL 

alpha 

Z3PAN 
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N»ITC  (6,906)  CPCPB 
WRITE  (6,916)  CPCHIT 
WRITE  (6,909) 

WRITE  (6,910)  CLIFT, CMOM 
GAMPRT«2.*GAMTE(i)*CPCPB 
WRITE  (6,908)  Z ( 1 ) , GAMPRT 
WRITE  (6,912) 

WRITE  (6,915)  (X ( I ) , Ta ILE , 1TE ) 

WRITE  (6,913) 

WRITE  (6,915)  (0(I),!«ILE,ITE) 

WRITE  (6,919) 

WRITE  (6,915)  (C ( I ) , I«ILE , ITE) 

900  FORMAT  ( 1 HI ) 

901  FORMAT  (30X,8A10) 

902  FORMAT  ( 1 H , /, 1 H , /, 1H  ,/,*  MACH  NUMBER  ■wEl3,5) 

901  FORMAT  (*  SIMILARITY  PARAMETER  (X)  ■*E13, 5) 

904  FORMAT  (A  THICKNESS  RATIO  baEIS.S) 

905  FORMAT  (a  AIRFOIL  ANGLE  OF  ATTACK  (RADIANS)  ■AE13.5) 

906  FORMAT  (a  CP  SCALING  FACTOR  (CP/CPBAR)  baE13.5) 

907  FORMAT  (a  WING  ASPECT  RATIO  **£13,5) 

908  FORMAT  (1H  ,/,lH  , /,21)(AAtRF0IL  SPANWISE  COORDINATE  baEH.5 
1 a SECTION  LIFT  COEFFICIENT  «*E13,5) 

909  FORMAT  (1H  ,/,lH  ,/,*  AIRFOIL  FORCE  COEFFICIENTSA) 

910  FORMAT  (1H  ,/,3XaLIFT  «*E1 3,5 , /, SX AMOMENT  ABOUT  (X«0)  ■*Et3,5) 

912  FORMAT  (1H  ,/,lH  , /, 3X*A IRFOIL  STREAMWISE  COORDINATE*) 

913  FORMAT  (1H  ,/,lH  ,Z,1X*AIRF0IL  PRESSURE  COEFFICIENTS,  UPPER  ■* 

914  FORMAT  (1H  , /, 3X*AIRF0IL  PRESSURE  COEFFICIENTS,  LOWER  ■*) 

915  FORMAT  (3X10E13.5) 

916  FORMAT  (a  CRITICAL  PRESSURE  COEFFICIENT  (SONIC)  a*El3.5) 

RETURN 

END 

FUNCTION  FUP(XX,ZZ) 

REAL  KCAP, MB , I W ING 

COMMON  /DELTA/  DX (40) , 0 Y (40 ) , DZ (20 ) , AX l (40 ) , AX2 (40 ) ,Bxl (40 ) , 

1 8X2(40), CX(40),AYl(40),AY2(40),AZ 1(20), A Z2 (20 ),X(40),Y(40), 

2 Z(20),FPU(S0O),FPL(B00),PHluB(800),TM,lMl, JM,JMl,KN,KMl, JW, 

3 JWP1,JWM1,ITE,TLE,KSPAN,DVBU1 , 0YBU2,DYBL1 ,DYBL2, ALPHA , DEL, MS 

4 GAM,KCAP,NDB,TITLE(8),DOU8,IWlNG,ZSPAN,NKPRT,KPRT(20i 
C airfoil  upper  surface  slope  DISTRIBUTION 

FUP»-4tA(xx.,5) 

C DOUB  IS  THE  DOUBLET  STRENGTH  DUE  TO  THICKNESS 
DOUB*l ,333333333* ZSP AN 
RETURN 
END 

SUBROUTINE  GAMFUN 
REAL  KCAP,ms,iwihg 

COMMON  /DELTA/  DX(40),DV(40),DZ(20),AX1(40),AX2(40),BX1(40), 

1 Bx2(40),CX(40),AY1(40),AY2(40),AZ1(20),AZ2(20),X(40),Y(40), 

2 Z(20),FPIJ(800),FPL(800),PHIUB(800),IM,IM1,JM,JM1,KM,KM1,JW, 

3 JWPI,JWM1,ITE,ILE,KSPAN,0YBIJ1,DYBU2,DYBL1,DYBL2, ALPHA, DEL» MS 

4 GAM,KCAP,NDB,  TITLE (8) , DOUB, IWING,ZSPAN,NKPRT,KPRT(20) 

COMMON  /Gamma/  GAMTE1(20),GAMtE(?0),PGFF,GAMFF(20) 

DO  10  1*1 , KSPAN 

GAMFF  ( I ) bGAMTE]  ( I ) TPGFFa  (GAMTE  (I)«GAMTE1(D) 

GAMTEKI)bGAMTE(I) 
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10  CONTINUE 
return 

END 

SUBROUTINE  INIT4L 

i ssiipss 

_ • ® 1 i INI 

10  CONTINUE  1 

£°  ^0  i«i,jmi 

*0  CONTINUE  1 

g°  JJ  I«1#KM1 

10 

"'MW  «*w,e,ttTt 

«0  CONTINUE 
J°  50  I«2,jmj 

JO  60  I«2,KM! 


6IREOIL  Bnnum.u  li 


■et  

..  KF 

▼o  CONTI Nue 
return 
eno 
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COMMON  /DELTA/  DV(40),DY(40),DZ(20),AXi  (40),  AX2(40),BXi  (40), 
a (a05  * **2(40), AZi  (28)  , AZ2(20),X(40),Y(40)  , 

f ZfiJ?)^PU(,025#FP,-(fl0n)»PHIuB^80MM»IMi»JM,JMJlKMlKMl,Jw, 

3 JWP1,JWMJ , ITE,ILE,KSPAN,DyBU1»DYBU2,DyBL1»DVBL2, ALPHA. DEL.  Mg. 

« «*m^CAP>NDB>TITLE(B).DOUB,IW1MG,Z3PAN,5kpRT,KPRT(20) 

COMMON  /COEFF/  A(40), 6(40), C(40),D(40),  PH 1(11500) 

COMMON  /GAMMA/  GAMTEl (20),GAMTE(20),PGFF,GAMFF(20) 

K9PAN1«K3PAN«1 

DO  10  K«1,KSPAN,K3PAN1 

MMl"ijM*(Kil  )E*n"1,,*tSAMfF(K>"GAMTEtK,)/(X(IMn"1,)*SAMTE(K,) 

IJK»MP*ITE*JM*JW 

PHI(IJK)»PHI(IJk)-PART 

L«MP*(ILE-2)*JM*JW 

LP«(ILE.2)*KM*K 

PHIIjB(LP)«PHI(L) 

LP«ITE*KM*X 

PHIUB(LP)«PHI(IJK)»2, OPART 
COMPUTE  CP  LOWER  (A)  AND  CP  UPPER  (§) 

DO  20  I«ILE,ITE 

L«MP*(Nl)*JM*Jw 

LP«(I-1)*KM*k 

!!5!"’!,*i!,<1C.n*CPNI{‘-4JM).PH|(L))^AX2(n*(PHl(L).PHICL.JM))) 

J JJ^CPHKjBfLP^KMl-PHTUBJLP)  )*AX2( I )* (PHIuB(LP)- 

1 PHIUB (LP»KM) ) ) 

20  CONTINUE 

PH1(IJK)«PHI(1JK)»PART 
WRITE  (6,901)  NITERG.K 

(B(I),I«1LE,ITE) 

NJTERG#  K 
(A(I),I«ILE,ITE) 


ITERATIONM5* 
TO  ITE)  ■ *) 

ITERATION*!?* 
TO  ITE)  •*) 


AND  K »*I3*  SCALED  PRESSURE  COPFF 


AND  K »*13*  SCALED  PRESSURE  COEFF 


WRITE  (6,902) 

WRITE  (6,905) 

WRITE  (6,902) 

10  CONTINUE 

901  FORMAT  ( 1H  ,/,*  AT 
1ICIENT,  UPPER  ( ILE 

902  FORMAT  (10E13.5) 

905  FORMAT  (1H  ,/,*  AT 

1ICIENT,  LOWER  ( ILE 
RETURN 
END 

SUBROUTINE  TRI  (I,K) 

REAL  KCAP,MS, IWING 

COMMON  /DELTA/  OX(40),Dv(«0),DZ(20),Axt(«0),AX2(iO),BXl(60), 

1 Bx2(40),CX(ttO),AYl(a0),AY2(«0),AZl(2O),AZ2(20),X(40),Y((lO). 

2 Z(2O),PPU(BOO),FPL(0OO),PMlUR(SOO),IM,lMl, JM, JM1 ,KN,KM| ,Jw, 

« iI^:i!!!,:IIf:Jtr:!JJp*;»0v®u,»0v®u2»0v8L1»t>vB‘-2»A*.PHA,DEL,MB, 
*.«®A!I,KCiPiN0*,T,TLe<®,,DOU,,»I“lNG**8>AN*NK^*T,FPRT(20) 

COMMON  /COEFF/  A{40), 6(40), C(40), 0(40), PH 1(11500) 

MP«IM*JM*(K«1) 

DO  10  KK*5,JH1 

J»JM1«KK*3 

P«A(J.l)/B(J) 

B(J-1)«B(J-1)*P*C(J) 

D( J»l )*D( J«1 )«P*D( J) 

10  CONTINUE 

MMMP*(I»1)*JM 
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00  20 
L«M+J 

20  CONTINUE 
RETURN 
End 


# 


-79- 


APPENDIX  B 


FORTRAN  LISTING  OF  TDUTRN 


A FORTRAN  listing  of  the  source  deck  for  the  TDUTRN 
program  is  presented  in  the  following  pages.  The  program, 
as  configured  here,  requires  161. 7aK  words  to  load  and  150. 0aK 
words  to  execute.  In  this  configuration  the  programs  fit  into 
small  core  of  the  CDC  7600. 
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2 

3 

« 

5 


C 

C 

c 

c 

c 

c 

e 

c 

c 


1 TAPF8, TAPER)  (INPUT*°UTRUT#TAPE5«INPUT#TAPEB»0UTPUT,TAPE7, 

1 0ME62l!siGI^pi,pT!?PuP?PLGTPHlIETBGtMTE'6AMFF'CRR'ERRnR» 
real  kcap.Js  * j,ppl,tphir»TPHIl»t^ill»tphibk,con 

2(20),PPj(800)fFPi!80fli  PHT.iJ!iiJ!,f2,!AZ2t20,»>((<,0>»Vt«0), 

NOOUfl.CPCPB.TITlE  8 ,Ie  OEL  A!!BL!*9MALLK»0Mt®» 

Z3PAN,KCAP,RPAR  >»"°»DEL, ALPHA, ITYPF,I0PT,XP,n.<PRT,KPRTC20). 

COMMON  /steady/  PHIsii?Sftft?A*!fi??2fP6FF,GAMFpC*0> 

i CX3U0),PHIuB8(*00) 1 * lSU0),Ax23U0),Bxi8(a0),BX23U0), 

COMMON  /INTEPP/  ZE(2SJ,N7E 

1 OMEGAE^MEBAP.rPSBPolNDUM^ikpiD  J«!,I2EiN,6AHFF»OHt6AM» 
a *prt, small*. i;,xp^?;J?^;;^ 

"AMELIST  /CONTRL/  mpE  t zc*  zc 

TO  ST^RT  PROCRlMi  iff  Any  iqu  tn  A A _ . _ 

’*«■•  unsteady  data  S.u,^,.r0‘i:,??e;noSeT:fE*?/'0"  * me 

NON  BE  WRITTEN  ON  TAPE7,  ILt  TAPE,«  NCw  UNSTEADY  DATA  WILL 
READ  STEADY  SOLUTION 

1 KSPANS#KCApJdEL  * ALPHA ^NDB^MS* BAM* iTESelLES, 

«CAD  (8)  (DUM,DUM,AXIS(I),AK2aJi{  *5yij?*MiDi!M,,>UM»DUM»z®PAN 
READ  (S)  OUM  *I3CIJ»AX2SCI)»BX1S(I), Bx28(I),CX8(I) * Ib | # IMS) 

READ  (S)  dum 
L«ITES*KMS 

«js  m 

L»IMS*JM3*WM3 

read  (8)  (PHIS(I). lal.Ll 

M00MiMwsr  trmling  coct  fhi 

MClw(lLEI«n*JMSTjws 

HC2«(ITES»n*JMS»jWS 

mci«(ites«2)#jms*,;ws 

LRl*(lLE8-n«RM8 

lR2«(ITES-n*KM8 

LP3«(ITE8-2J*rms 

00  1 J , kspans 

MPbh«(k.1) 

L«MPtmC1 

LRwLPItk 

L*MPfMC2  ^Ml8  ^E)^PH|uB8  (LP)  ) 

LP*LP2*K 
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v_ 


c 


c 

c 


c 


c 


PHI3(L)».5*(PHIS(LWPHluB3(LP)) 

L»Hp+MC3 


LP»LP3*K 

PHI3(L)»,5*(PHI3(L)TPHIUB3(LP)) 
1 CONTINUE 

SKbSQRT(KCAP) 


CPCPBaDEL**, 6666666667/ ( (l,+GAM)*Mft**2)**( 3333333333 
R p A R a } f / ( (J (+GAM)*MS**2*DEL)**, 33333333 3 3 
READ  (5, PH)  (TITLE ( I ) , 1*1 ,6) 

READ  (5,CONTRL) 

IP  (ITAPE.EQ.O)  GO  TO  10 
READ  DATA  PROM  RESTART  TAPE 


READ  (9)  (PPU{I),FPL<n,PHlUBCn,I«l,l) 

READ  (9)  (GAmTE(I),GAMFF(I),Xb1,K3PAN) 

L»IM*JM*KM 

READ  (9)  (PHI(U(I«i,Ll 
00  2 IB1,KSPAN 
GAHTE1 ( I JbGAMTE ( I ) 

2 CONTINUE 
IK»0 

READ  (5, IN) 

WRITE  (6,900) 

WRITE  (6,901)  NITERG 

niterg»o 

WRITE  (6,913)  KCAP,CPCPB 

THE  IK  OPTION  IS  U3FO  TO  BOOTSTRAP  TO  DIPFERENT  REOUCED  FREQUENCIES 
AND/OR  MOOES  OF  OSCILLATION 
IF  (IK.EQ.O)  GO  TO  15 

0MEGbSMALLK*M6*«2/( ( 1 , *G AM )*DEL*M6* *2) **.6666666667 
CALL  INITAL 

call  farflo 

GO  TO  15 

START  PROBLEM  PROM  SCRATCH 
10  CONTINUE 
READ  (5, IN) 

NITERG»0 

NDOUBaO 


OMEG»8MALLK*M6**2/( ( 1 ,*G AM) *DEL*MS« *?)**, 6666666667 

00  3 I*1,K3PAN 

CAMTEICI)bGAMPP(I) 

GAMTE(I)»GAMPPCI) 

S CONTINUE 
IM1bIM«1 


JmjbJm»i 

KMJbKM-1 


JWPlajWTl 

JWMlaJWl 

INITIALIZE  FINITE  DIFFERENCE  COEFFICIENTS  AND  FARTIELO 
CALL  INITAL 
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CALL  FARPlD 

INITIALIZE  GUESS  FOR  SUBSONIC  CASE  (INTERIOR  ONLY) 
ASSUMED  SYMMETRY  IN  Y U 

Jw2«2*Jw 

CON«GAMFF(l )/6, 2631853 
DO  20  K«1,KM1 
MP«IM*JM*(K-1) 

ZP«Z(K)*Z8PAN 
ZM«Z(K)-Z3PAN 
Z2«Z(K)**2 
DO  30  I«2, IM1 
M«MP+(I-1 )*JM 

PHI(MTjW)«CMPLY(fl,,Ot) 

X2«X(I)**2 
do  ao  JbJwpi,jmi 
L«mtj 

LL»m+Jwz*j 
Y2«Y ( J) **? 

IF  (XdJ.LT.l.)  GO  TO  «1 

PHI(L)«CON*(ATAN(ZP/Y(J))-ATAN(7M/Y(J))) 

CO  TO  02 
01  CONTINUE 

R»3QRT(X24KCAP*CyZ»Z2)) 

„ !S!S"W!P‘ 

CPHI«CAflS(PMI(LJ) 

IF  (CPHI.GT.l.)  PHICL)“PHI(L)/CPHI 
PHI(LL)«-PHI(L) 

40  CONTINUE 
SO  CONTINUE 
20  CONTINUE 
L*ITE*KM 

ERR«CMPLX(0,,0,) 
do  a i«!,l 

PHIUB(I)«ERR 
a CONTINUE 

N»(ILE-2)*JMfJW 
KK«(ILE-1)*KM 
DO  05  K«1#K3PAN 
L«MtIM*JM*(K-J) 

PHIUB(KK*K)«PHI(L) 
as  CONTINUE 
IS  CONTINUE 

OMEG2I«CMPLX(0,,2,*OMEG) 

WRITE  (6# IN) 

WRITE  (6*900) 

IF  (ITAPE.EQ.O)  WRITE  (6,913)  KCAP,CPCPB 
KCRD«1 

RI-CYCLE  POINT  FOR  GRID  ITERATION 
SO  CONTINUE 

ERROR«CMPLX(0,,0,) 

NIT«NITERG 

NITCRCBNITERCTI 

IF  (MOD(NITERC,NPRINT).EQ.O)  CALL  PRINT (NIT) 

IF  (MOD(NITERC,NGFF) ,NEf 0)  GO  TO  51 
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51 


HITERG,  GAMTE  n ) rautrio 

cm,,, 


PLANES  (Z  DIRECTION) 


102 


•1) 


call  gamfun 
C*LL  EArfL0 
*"ITC  (6,910) 

COMTIMIJE 

C6<0B«6»DOU8) 

K8p1*INCR 
I HjMajH* JM 
JM'5JW8«IM3*JH3 
REGIN  LOOP  ON  THE 
00  100  Kbj , kmj 
K8pK3MNCR 

00  102  1*2, JMI 

VfljMKCAP 
CONTINUE 
HP«lMjM.(K„n 

CHECK  FOR  AIRFOIL 
I^OIL«0 

,rn,L" 

'"so  tsr.3.^”**-  ’U,,c  '*  OI'fCTION) 
ia.mincR 
CHECK  FOR  AIRFOIL 

iklag-o 

..uUS,M^5*(l5-n*jM8 

# 5«TW,S  column  of  PHI 
00  201 
L»H9J 

•'HIogtjjbpmjj,. 

CONTINUE 
J8«2*INCR 

8CQIN  LOOP  ON  COLUMN  fv 

»o  m ,v  t,,,«noN) 

...  f««J8tINCR 

e\eftii:„e,“  ,noic"  «» «... 

L8RPL89JM8 

L8L«LS-JM8 

LSLL«LSL-JM3 

LR»Lpjh 

li«l»jh 

Ul*LL»JH 
;;  n.Eo.?, 

LBpl.1 
tA»LPl 
i-^PLPlHJM 

LRKpl-Ihjm 
^ fK.EG.n 


IPLAG*1 


201 


ULpIL 


L8K«lF 
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nnn  n n n n n n tin  t» 


* CALCULATE  V AND  PHIXX  FROM  STEADY  SOLUTION 


VVSKCAP.AXtS(IS)*(PHla(LaR).PMIS(LS))-AX2S(IS)*(PMI3(L3). 

I PMlStLSL)) 

WS«VV 

IE  (VV.LT.O.)  60  TO  SOI 
ELLIPTIC 

OMEGA ( J)sQMEGAE 

PMIxXbBX13(IS)*(PHIS(L3P)-PMIS(L3))-BX23(IS)*(PHI3(L3).PHI3(L3L)) 
60  TO  302 

501  CONTINUE 
OMEGA (J) IOMEGA P 

IF  (V(J).GT.O.)  GO  TO  SOS 
HYPERBOLIC 

OMEGA (J)bQMEGAH 

VViKCAP.CXS  (IS-1)*  (PHIS (L3)-PHIS(L3D) -CXS( 18-2)* (PHIS (LSL)- 
1 PMIS(LSLL)) 

SOS  CONTINUE 
PARABOLIC 

PHIXX»BXI S(IS-l)* (PHIS (LS) -PHIS (L8L))-BX2S(IS-I )* (PHIS (L8L)- 
l PMIS(LSLL)) 

502  CONTINUE 
V(J)aVV8 
TPHIRsPHI(LR) 

TPHIL«PHI(LL) 

TPHILL"PHI(LLL) 

TPMIBKbPHKLBH) 

IE  (lEOlL.EQ.O.OR.J.NE.JM)  GO  TO  SOS 

IE  (I.EQ.ILE-I)  PHI(LR)*,5*(PMIUB((ILE-1)*KM*K)*PHI(LR)) 

IE  (I, EG. ITEM)  PHI(LL)*,5*(PMIUB((m-I)*KM*K)*PHI(LL)) 

IE  (I, EG. ITEM)  PHI(LLL)*.5*(PHIUB((ITE-2)*KM*k)*PHI(lLL)J 
IE  (I,E0,ITE*2)  PHI(LLL)*.5*(PHIU8( (ITE-1 )*KM4K ) 4PHI (LLL ) ) 

SOS  CONTINUE 

IE  (ILC,LE,I,AND,I,LC,ITE,AND,J,EQ,JM,AND,K,EQ.K3PAN*1) 

1 PHZ(LBK)*,B*(PHIUB((I-1)*KM*K3PAN)*PHI(LBK)) 

BET  UP  TPIDIAGONAL  MATRIX  TO  SOLVE  FOR  PMI(I,J,K) 

A * PHI(I,J*I,X)  ♦ B * PHId.J.K)  ♦ C * PHI (IfJ-lfK)  ■ D 
IE  (IPLAG.fG.I.ANO.J.EQ.jNPl)  GO  TO  SSO 
IE  (IPLAG.EQ.I.ANO.J.EQ.Jw)  GO  TO  SSO 
IE  (lELAG.EQ.I.AND'J'EQ.JMM!)  GO  TO  550 
PART«CHPLX(0,,0.) 

IE  (I.LE.ITE.OR.IEOIL.CO.O)  GO  TO  305 
KUTTA  CONDITION 

S!G!«(X(I)-t,)*(GAMFF(K).GAMTE(K))/(X(IM!).i#)4GAMTE(K) 

IE  (lOPT.EO.I)  S1GI*8IGI*CEXP (CMPLX (0, * -SMALLH* (X ( I ) -1 , ) ) ) 

IE  (J.EO.JHMl)  PART«,5*AY) (J)*SIGI 
IE  (J.EQ.JM  PARTS, S*(AV1(J).AV2(J))*SIGI 
IE  ( J,EQ« JNP1 ) PARTS., 5«Ay2(J)*SiGl 
SOS  continue 

IE  (VVS.LT.O.)  GO  TO  320 


* ELLIPTIC  DIFFERENCING 


A(J)sAY1(J) 
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ooo  r*  r»  n nn 


B ( J J * 


1 


iiTT^?n2?15i?:2JiJir*AvlCJi+AV2<JWC0MC02I^MIX,<)*C4,<2C1)‘ 


i 


C(J)«AY2(J) 
D(J)« 


if  CJ.co.jMt)  go  to  M2 

CO  TO  390 
BOTTOM  BOUNDARY 

511  CONTINUE 
0(J)»0(J)-AY2(J1*PHKLB) 
GO  TO  190 

TOP  BOUNDARY 

512  CONTINUE 
D(J)»D(J).AYKJ)*PHI<lai 
GO  TO  J90 


HYPERBOLIC  AND  PARABOLIC  DIFFERENCING 


320  CONTINUE  " »•••••••******, 

*( J1«AY1 CJ) 

1^4Z2IK)^*^  ^I"I5"^^1IJ)"*T2(J) * COME G2 It PH! XX ) *CX(I«l)"AZl  (K ) • 
C(J)«AY2(J) 

0(J)*VV*(BX1 (I-1)aPHI(LL5*Bx2(I-1 > * (PHI CLL >»PHI (LLL) ) )-(OMEG?I* 

IF  (J.E0.2)  GO  TO  321 
IF  (J.E0.JM1)  GO  TO  522 
GO  TO  390 
BOTTOM  BOUNOARY 

321  CONTINUE 
D(J1«D(J)-AY2(J)*PHI(LH) 

GO  TO  390 

TOP  BOUNOARY 

322  CONTINUE 
D(J)«D(J)-AY1 (J)aPHI(LA) 

GO  TO  390 

************  *IFFOIL  UPPER  SURFACE  BOUNDARY  CONDITION********** 

330  CONTINUE  ***************************************** 

IF  (VVS.LT.O.)  GO  TO  331 
ELLIPTIC 

*( J)»DYBU1 

I8C  AZicij-Izicij15  f8x2<15  5*0vBui^<OM<r62I^,'HIX,<5*CAX2CI5-Axl  <15)3- 

C(J)*0, 

D(J)«0YBU2*FPU(N).w*(BXl  (I)*PHI (LR)*BX2( I )*PHI (LLJ) ♦IOMEG2I* 

GO  TO  190 

HYPERBOLIC  AND  PARABOLIC 

331  CONTINUE 
A(J)«DYBU1 
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A 


8(J)«VV*Bxl(I-n-PVBUl-<0MEG2I*PHIXX)*CXtI*n-AZl(K).AZ2(K) 

CtJ)«0, 

0(  J)«DYBU2*FPU(N)  * VV*(BXl  ( !•!  )*PHI  (LL)*BX2(I»D*(PHI  CLD* 

1 PHICLLU)  )-(OMtG2I*PMIXX)*(CXCI»n*PHl(LL)-CXCI»2)*(PHI  CLL >- 

2 PMiaLl)>)-fAZi<K)*FHI(LF)*AZ2<K)*PHICLB»<)  j 
GO  TO  390 

**  ***  AIRFOIL  LOWER  SURFACE  BOUNDARY  CONDITION 


HO  CONTINUE 

IF  (VVS.LT.O,)  GO  TO  351 
C ELLIPTIC 

B(  j?l!<DYHll*VV*(BXl<IHBX2<imtOMEG2I  + PHlXX>*(AX2(I)-AXl(n>)' 

1 AZl(K)«AZ20O 

D(j)«»DY8l2*FPl(N)»VV*(BXl (I)*PHI (LR)*BX?(I)*PHI (LLI )*(OMEG2I* 

1 PHIXX)*(AX1  (n*PMl(LR)-AX2tI)*PHT  ClL)>-CAZl(K)*PHI(LP>* 

2 A22(B)*PHI(LBK)) 

GO  TO  390 

C HYPERBOLIC  AND  PARABOLIC 
HI  CONTINUE 
A(j)-0. 

B(J)«VV*8XlCI-n-OyBLl-(OHEG2IyPHIXX)*CX(I-l)*AZl(K)-AZ2(K) 

C( J)«DYBL1 

D(  J)b»DYBL2*FPL(N)aW*(BX1  € I-i ) «PMI  (LL)+BX2( I-D*  (PHI  (LL)- 

1 PHI(LLL)i)-(0MEG?I*PHIXX)*(CX(I-n*PHI(LL)-CX(I»2)*(PHI(LL)- 

2 PHl(LLUn-(AZl(K)*PHI{LF)TAZ2CO*PHl(LBK)) 

GO  TO  390 

C BODY  BOUNDARY  J«JW 
340  CONTINUE 
A<J)«0, 

B(J)»CMPLX(1,,0,J 

C(J)«0, 

D(J)«PHI(L) 

140  CONTINUE 

PHI(LR)«TPHIR 
PHI(LL)»TPHIL 
PHI (LLL)"TPHILL 
PHICLBK)«TPHIBK 
IF  (IOPT.EQ.O)  GO  TO  300 
IF  (IFLAG.EO.t.AND.J.EO.JN)  GO  TO  300 
8{J)nB(J)4SNALLK*0HEG 
300  CONTINUE 

C TRXDIAGONAL  matrix  IS  SET  non  solve  for  COLUMN  of  PHI 
CALL  TRI (I»X) 

C RELAX  phi,  find  ERROR  and  move  to  next  column 
DO  345  JB?,JM1 
L»M*J 

ERRaQMEGA ( J)*(PHX (l)»PHIOG( J)  ) 

PHI(L)«PHIOG(J)AtRR 

IF  (CABS(ERR) ,LT,CABS(ERROR) ) GO  TO  345 

errorrerr 

LERROR«L 
341  CONTINUE 
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1 


IP  (IPLAG.NE.l)  60  TO  200 

!?IV5C!1*?;L(L:1)*0YCJW,,,^HI{^2VPHICLM  /OY  j!pn 

m 

too  continue 

: PRINT  OUT  EPRCIR  APTPR  EACH  GRID  SWEEP 

WRITE  (6,9051  NITERG, ERROR, LERROR 
IP  (CABSCERRI.LT, 10,1  GO  TO  101 
WRITE  (6,912)  1 

STOP 

101  CONTINUE 
IDOUBaO 

IP  (CABS(ERROR),LE,tPSGRO(KGRD1)  GO  TO  900 
IP  (NITERG, E9.NGRID)  GO  TP  fllO 
IP  (MOOtNITERG.NDUHPI.Efl.O)  GO  TO  910 
GO  TO  50 
•00  CONTINUE 
KGRDaKGRDM 
IDOUBal 
GO  TO  910 
•01  CONTINUE 
CALL  GAMPuN 

CALLCpPRiNT0)  *lTERe,GAMTE(1)»CAM*Pn),GAHTt(KSPAN),GAHPF(KSPAN) 
WRITE  (6,900) 

WRITE  (6,906)  NITERG 
CALL  DOUBLE 

I^JB,JWiKM.ILC.ITE,KSPAN 

WRITE  (6,902) 

WRITE  (6,903)  ( X (I ) , Ial , IN) 

WRITE  (6,909) 

WRITE  (6,903)  (Y(I),Ia}t JR) 

WRITE  (6,915) 

WRITE  (6,903)  (Z(I ) , Ial ,KH) 

GO  TO  50 
•10  CONTINUE 

(^)  N I TERG,  I M , J H J , JM,  JR  J , KM  , K H J , JW  , JWP 1 , JwMl  . f TP  fi  r 

1yO5JSAi';?M?5:?!?A!:l:,!:?Y?H?;?TBV2»?!!,‘'*»D^8L2,ND0uB.xp 


(Z(I),DZ(I),AZl(I),AZ2(X),Ial,KR) 

(PPU(I),PPL(X),PHIUB(I),IaJ,L) 

(GAHTE(X),GARPPU),X«1,KSPAN) 


write  (7)  (X(I),0X(I),AX1(I),I 
WRITE  (7)  (V(I),DV(I),Avi(I),I 
WRITE  (7)  (Z(I),DZ(I),AZ1(I),I 
L«ITE*km 

WRITE  (7)  (PPU(I),PPL(X),PHIU! 
WRITE  (7)  (GARTE(I),GARPP(I),] 
L*IN* JM*KR 

write  (7)  (PHI (I ) , Ial , L) 

END  PILE  7 

WRITE  (6,907)  NITERG 

CALL  PRINT(NITERG) 

IP  (WORO.GT.KEPS)  GO  TO  920 
IP  (NITERG, EG, NORXD)  GO  TO  930 
IP  (XDOUB.EO.l)  GO  TO  901 
GO  TO  50 
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420  CONTINUE, 

WRITE  (S,90«)  NITERS 
GO  TO  4*>0 

4130  CONTINUE 

WRITE  (B,9ft9)  NITERG 

R01  FORMAT  (1H  ,/,*  CASE  IS  BEING  RESTARTED  AT  ITERATION*!*) 

902  FORMAT  ( 1 H ,/,*  X(I).I«1,IH*) 

903  FORMAT  (10F13.5) 

905  FORMAT  (!m  A^1 I Tf RAT ION* 15*  THE  MAXIMUM  ERROR  ■•2C13.S*  AND  0 

909^  FORMAT^  (IN  (!  / # **THE  NUMBER  OF  NODES  IS  BEING  DOUBLED  AT  ITERATION* 

907%ORMAT  (IN  ,/,*  TAPf  HAS  BEEN  DUMPED  AT  ITiRAT  ION*I5) 

90S  FORMAT  (1M  ,/,«  SOLUTION  HAS  CONVERGED  TO  DESIRED  ACCURACY  AT  ITE’ 

ROr'fORmAT1 (1H  ,/,*  MAXIMUM  NUMBER  OF  ITERATIONS  HAS  BEEN  REACHED,  CAS 
IE  IS  BEING  TERMINATED  AT  ITERATION#I5) 

910  FORMAT  ( 1 H ,/,*  UPDATE  GAMpF  AND  FARFIELD  AT  ITERATI0N*!5,4X*GAMT 
,r,o  I-Jnt.i  4x*  GAMFF(l)  **201.5, /,44X*  GAMTE (KRP AN)  **2El3.5* 


¥0^  r UR  n • T v I n § f 9 " un 

IE  IS  BEING  TERMINATED  at  ITERATION#!*) 

910  FORMAT  (1H  ,/,#  UPDATE  GAMpF  AND  FARFIELD 
1E(1)  «*?tl3.5,4X*  GAMFF(l)  **201, 5, /.<*«*• 

a -ASCII  « 3 


2GAMFF(KSPAN)  **2E!l,5) 

911  FORMAT  (BA l 0 ) 

912  FORMAT  (1H  ,/,#  SOLUTION  IB  DIVERGING, 

1ATEO#)  . 


9U  FORMAT  (1H  ,/,#  SIMILARITY  PARAMETER  (*)  »*E13.5,/. 
1R  (CP/CPBAR)  «*E13,5) 

9)4  FORMAT  (1H  ,/,#  IM  **I  4*  JM  **I4*  JW  »*I4*  KM  ■•!«< 
1 # ITE  ■*IR*  K9MAN  #*14) 

915  FORMAT  (1H  ,/,*  2(I),!«1.KM#) 

450  CONTINUE 

CALL  FPRINT 
END 

SUBROUTINE  DOUBLE 

COMPLEX  B,  D,  PHI  IJB,  PHI,  GAMTE  1 .GAMTE,  GAMPF,  FPU,  FPL 

real  kc ap , mb 

COMMON  /DELTA/  DX(«0), DY(40), 02(40), AX1 (40), A*2C«0 

1 BX2(40),CX(40),  AVI  (40),AY2(40)  , AIH20)  ,A22(20)»* 

2 2(20) , FPU (BOO) ,FPl (BOO), PHI U§ (BOO) , IM, IMJ ,JH,JM1 
1 JHPl,jMMi,ITf.lLE,K8PAN,DYBlll,DYBU2,0YBLl»0YBL2, 
« NDQUB , CPCPB, TITLE (8). MB, DEL, ALPHA, ITYPE,IOPT,MP, 
9 ESPAN.KCAP.RPAR 

common  /COEFF/  A(40),B(40),C(40),D(40),PHI(11500) 
COMMON  /GAMMA/  GAMTE 1 (2ft ), GAMTE (?0 ), PGFF , GAHFF (20 ) 


, THE  PROBLEM  IS  BEING  TERMIN 
(K)  P*E13,5, /,#  SCALING  FACTO 
■* ! 4#  KM  I 4#  ILC  »#I4 


COMMON  /GA 

RETURN 

ENO 


COMPLEX1 B^D,PHltpHluB, FPU, FPL, GAMTE 1, GAMTE, GAHFF, PI ,PlO, PARTI , 

PARTI 0,OMK, WING, AMUK, WAKE  IN, Gi , G2, GAMJEI , CON4, CON5 

-k  i . A ■ U A 


5 7SPAN,KCAP,RP*R 

Ff  !iS  * 

,u”2.NiS4J«ip^  (*S,UME0  S¥MMET,,¥  in  »> 

C0Nl«l#/6,283iaS3 

CON2»KCAP*C0N1 

C0N3«0MEG/KCAP 

AMU*SQRT(0MEG*(C0n3+SMALLK11 

AMUK*CMPLX(0,#AMU/8K) 

8ETA2*1 

0mk«Cmplx(O.,-CON3) 

JW2»2* Jw 

IM JM»  JM# JM 

CALCULATE  PART  OF  WING  INTEGRAL 

on  10  i«ile,ite 

*L«(IM)*JM^JW 

MU«(I-1)*KM 

C0Na«CEXP(0MK*X(I)5 

PART1«CMPLX(0,,04) 

Do  20  K*1 , K3PAN 
L«IMJM*(K-1)*ML 
LP«MUfK 

P1»PHIUB(LP)*PHI(L) 

IP  fK.EG.n  GO  TO  21 

co"i”£"m,5*tP,*P,0,*OM'<*n 

P10«P1 
CONTINUE 

?"n:tTRE)‘"0T1;55;r*("p‘N*z(Ksp^>>> 

eoNT^^lueS*'S*tP,RT1*P*,,T'0,*0*tl*l, 

MUTt'oMiJrf’  “ING,*S*,’*f'Tl.»«ILE) 

CONTINUE 
INTEGRATE  GAMTE 

GAMTEI«CHPLX(0,,0,) 

00  IS  K»2#K3PAN 

E.Eti2)EI'G‘HT'I*,S*G‘MTE(l<SP*N,*t“,,*N-*0<3P*N)} 
MP«IMJM*KM1 
Z2>Z(km)**p 
oo  30 

M»HP*(I«1 j#jm 
X2«X (I)**p 
XO«X(I)-i, 

X0?sX0**2 

PHI(M4Jw)«CHPLX(0,,0.) 


21 

20 


11 


10 
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L«M*J 

LL"M*jw2-j 

V2«V(J)**2 

R«KCAP#(V2tZ2) 

BRa30RT(X02*R) 

IRbSQRT(R) 

Tla(M8#B«*X0)/BETA2 

RHaSR*RPAR 

U«Tl/RH 

CALL  WAKF  (U#3mALLK,RH#wAKEIN) 

GlaKCAP*MB*CEXP(CHPLX(0,«»3HALLK*Tl ) ) / ( BR* (BR»MR*XO ) ) 

62®WAKt I N/R 

8R«3QRT(X2*KCAP*CV24Z2) ) 

PHI (L)«COMR«V( J)*{G1*S2)4C0n5*V ( JJ* (| .♦AMUK*BR)*CEXP ( 
1 BRooJ 
PHKLL)a-PMI(L) 

si  continue 
so  continue 
c x»x(n 

X2«X  ( 1 )**2 
X0aX(l)-t, 

X02aX0**2 

CON<iaeON|*CEXP(CMPLX<0.,.SMALLK*XO))«GAMTEI 
CON5aCON2*alNG*CEXP(CMPLX(0,,CONS#X(i ))) 

00  60  Ka| , KM1 

Ha  I M JM* ( X » 1 ) 

Z2«Z (K)**2 

PHI(H*JWJaCMPLX(Of  ,0.) 

DO  61  JaJwPl'JM 
L»M4j 

LL«H4JW2*J 
V2avCJ)**2 
RaKCAP*( V2*22. 

BRasQRT (X02PR) 

IRaSORT(R) 

Tla(M8*BR«X0)/BETA2 

RHasR*RPAR 

UaTl/RM 

CALL  WAKE  (U,3HALLH*RH,wAKEIN) 

JJjJJJJjJJJ^XPfCHPLXtO^.SHALLKwTm/tBR.IBR.HBaxOn 

8Rb30RT(X2pKCAP*CY2pZ2)) 

JpCON5*V(J)*(1  .♦AMUR*BR)#CEXP(* 

RMI(LL)a*PNICL) 

61  CONTINUE 
60  CONTINUE 
C XBX(IM) 

XJaXMl *JM 
X2«X ( XH)««2 
XOax ( I M) ■ 1 , 

X02ax0**2 

CON«acON!*CEXP(CMPLX(0.,*3HALLX*XO)) 
CON5aCON2*alNG*CPXP(CNRLX(0.,CON5*X(IM))J 
DO  43  KalfKM 


■AMUK*BR) / 


AHUK*BR) / 
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1 


c 


)*IJ 

PHI (M^JW)»CMPLX(0,,0,) 

Z2«ZCK)**2 
00  70  I * 1 , *^2E 

A(n*zEm 

70  CONTINUE 
NEND«N*E 
IFLIP-0 

00  aa  J*  ,i  w P i (jm 

M(«Ji;uijp;2'(,Ifap-E0-n  00  70  71 

‘(Dim) 

72  CONTINUE 

neno»kspan 

IELIP«t 

71  CONTINUE 
L«m*J 

LLPM4JW2.J 

Y2«Y(J)**2 

PART1«CHPLX(0.,0.) 

1 2 E * 2 

00  a?  KK«1,NEN0 
a7  CONTINUE 

IP  CA(KK),LEaZ(IZE))  GO  TO  48 
IZE«IZEM 
GO  TO  az 
48  CONTINUE 


R«KCAP*CY2*CZ(kj-a(kk))**2) 

BR«3QRT(X02tR) 

9R»SQRT(R) 

T1«(H8*BR.X0)/BETA2 
RH«3R*RPAR 
U«T1 /RH 

C*LU  WAKE  CU*SHALLK,RM#WAkEIN) 



IP  (KK.EO.l)  go  to  a 6 

PARTI«PARTl4,«j*{pltpln),(AUKj.A(KK- 

CONTINUE 

PlO»Pl 

CONTINUE 

?H{tU)«pIl(UN5*V<J,‘<,'’*HUK*(i,,*CE*,,<**H,IK*l,R,/S''**1 
a a continue 
aj  continue 

VRYCl)  AND  Y«Y(.JH) 

Y2«Y ( J) **2 
00  S3  K«1,KM! 

MP«IMJM*(K*1) 

Z2«Z(K)**2 


46 


as 
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•MB*XO)) 


00  59 

)*jm 

L«M  + J 

LL«M*JW2*,J 
X2«x(I)**2 
X0«X(I)-1 , 

X02«X0**2 

CO^a«CnNl*CEXP(CHPLX(0..*9MALLK*yft1wr*MTc  t 

52KNc5:^:;rcE*,,<c’,BU(c"cnNj-*ti'» 

BR«SQRT(X02*R) 

3R>9QRT(R) 

T1«(M8*8R.xO>/B£TA? 
rh«sr*rpar 

U»T  1 /RH 

CALL  WAKE  (U, SMALL*, RH,wAKEIN) 

8R«SQRT(X2*KCAP*(v2+Z?)) 
i»HI,uK0,..y(J).(cu«)tC0N5.r(J).(I,ttMUK,|R),ct<p(>iiiuK4M)/ 

PHI (LL)«-PHI(L) 

59  CONTINUE 
51  CONTINUE 
return 

END 

SUBROUTINE  EPRINT 

2 2(20),PPu{«0n),EPL(800),PMIuJ  sS5i  l?2,:Jf2**0,t:,,Ca0,*vC«°>. 

CPOEL-CPCPB/OEL 
WRITE  (6,900) 

WRITE  (6,901)  (TITLE (!),!a!,0) 

WRITE  (6,902)  MS 
WRITE  (6,903)  KCAP 
WRITE  (6,909)  OEL 
WRITE  (6,905)  ALPHA 
WRITE  (6,906)  SMALL* 

WRITE  (6,907)  OMEG 
WRITE  (6,900)  xp 
WRITE  (6,909)  Z9PAN 
WRITE  (6,910)  CPCPB 
WRITE  (6,911) 

i6l9,2>  (x<IM«ILE»ITE) 

CLIFT«CMPLX(0#,0,) 

CMOM«CMPLX(0,,0#) 


-93- 


0?  10  K«1,K9PAN 

PART*. 5*C(X (ITEM )-l,)*(CAMFFCK).GAMTECK))/(X(lMI).i.)*GAMTE(K)) 
IF  (IOPT.EQ.l)  PART*PART*CEXP(CMPIX(0.,»SHALLK*(X(ITE+1)*1,))) 

IJK*MP*ITF.*JMfJw 
PHI CIJK)«PHI (IJX)-PART 
L«MP4CLF-2)*JM*JH 
LP*  C ILE»2) *KMfK 
PHIUB(LP)«PH! (L) 

LP*ITE*K*HK 

PHIuB(LP1«PHI(!JK)42,*PaRt 
00  20  I*ILE» ITE 
M«MPf(I.l)*jM 
L«H4jW 

LP*(I-1)*KM*K 

B(I)*»2.*(AX1(I)*(PHX(L+JH)*FhX(L))+AX2(I)*(PHI(L)*PHI(L«JM)))* 

1 CPOEL 

D(I)*-2,*(Axi CI)*CPHIUB(LP4KH)*PHIUBU.P))*AX?(I)*(PHIUB(LP)» 

1 PHIUB(LP-KM)))*CP0EL 
IF  (inPT.EO.O)  GO  TO  2« 

C1«CHPLX(0.,2,*9MALLK)*CPOEL 

B(I)*B(I)-C1*PHI(L) 

D(I)*DCI)-Cl*PHIUB(LP) 

2«  CONTINUE 

IF  (R.GT.l)  GO  TO  21 
Bl(I)*B(I) 

Ot(I)*D(I) 

21  CONTINUE 
C1«B(I)-D(I) 

C2«Cl*(X(I)-XP) 

IF  CI.GT.ILE)  GO  TO  22 
CL«C1*XCILE) 

CM* ,5*C2*X  ( HE ) 

GO  TO  25 

22  CONTINUE 

CL«CL+.5*CClfC10)*OX(I-l) 

CH«CM4,5*CC2+C20)*OX(I-1) 

25  CONTINUE 
C10*C1 
C20*C2 
20  CONTINUE 

PHICIJK)«PHI(IJk) »PART 
IF  CK.EQ.l)  GO  TO  11 
CLlFT«CllFTf,5*(CLfCl.0)*DZCK-l) 

CMOH«CMOH+,5*(CH*CMf))*OZCK»l) 

11  CONTINUE 
CLO«CL 
CHO«CM 

00  12  N*1 , NKPRT 
IF  (KPRT(N).NE.K)  GO  TO  12 
GAMPRT«2,*GAHTF(K)*CP0EL 
WRITE  (6,913)  Z(K),GAMPRT 
WRITE  C6,9l«) 

WRITE  (6,915)  (0(I),I«ILE,ITE) 

WRITE  (6,916) 
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WRITE  (6.915)  (B(I),I«ILE,ITE) 

GO  TO  10 
12  CONTINUE 
10  CONTINUE 

WRITE  (6,909) 

WRITE  (6,901)  (TITLE  ( I) , I»1  #8) 

WRITE  (6#  9021  *8 
WRITE  (6,903)  KCAP 
WRITE  (6,90a)  OFL 
WRITE  (6,905)  ALPHA 
WRITE  (6,906)  SMALL* 

WRITE  (6,907)  OMF.G 
WRITE  (6,908)  XP 
WRITE  (6,909)  7SPAN 
WRITE  (6,910)  CPCPB 
GAmprt«2,*GAMTE(1) *CP0EL 
GO  TO  (30, 35,40) , ITVPE 
30  CONTINUE 

WRITE  (6,917) 
write  (6,9i8)  clift, chop 
white  (6,913)  Zm.GAMPHT 
WRITE  (6,919) 

GO  TO  45 
35  CONTINUE 

WRITE  (6,920) 
writt  (6,918)  clift, cmoh 
write  (6,913)  Z(1),GAMPRT 
white  (6,921) 

GO  TO  45 
40  CONTINUE 

WRITE  (6,922) 
write  (6,918)  clift, CMOH 
WHITE  (6,913)  Z(1).GAMPRT 
WHITE  (6,923) 

45  CONTINUE 

WRITE  (6.911) 

WRITE  (6,912)  (X(I),I«ILE.ITE) 

WRITE  (6,914) 

WRITE  (6,915)  (01(T),I«ILF.,ITE) 

WRITE  (6,916) 

WRITE  (6.915)  (B1(I),I»ILE,ITE) 

900  FORMAT  ( 1 H 1 ) 

Ofll  FORMAT  (30X,8A10) 

902  FORMAT  (1m  ,/,lH  ,/,*  MACH  NUMBER  ■*E13,5) 

903  FORMAT  (*  SIMILARITY  PARAMETER  »*E13,5) 

904  FORMAT  (*  THTCKNFSS  RATIO  ■*El3,5) 

905  FORMAT  (*  AIRFOIL  ANGLE  OF  ATTACK  (RADIANS)  «*C13.5) 

906  FORMAT  (*  REDUCED  FREQUENCY  (BASED  ON  CHORD)  ■*C13,5) 

907  FORMAT  (*  SCALED  FREQUENCY  (OMEGA)  »*E13,5) 

90S  FORMAT  (*  PITCH  AXIS  (XP)  «*Ei3,5) 

909  FORMAT  (*  WING  ASPECT  RATIO  «*E13,5) 

910  FORMAT  (*  CP  SCALING  FACTOR  (CP/CPBAR)  «*E13,5) 

911  FORMAT  (1H  ,/,lH  ,/,3X*AIRFOlL  STREAMwISE  COORDINATE*) 

912  FORMAT  (3xei3.5,13XEn,S,l3xEl3.5,l3XF13,5,l3XEt3.5) 

913  FORMAT  (1H  , / , 1 H , /, 15X*AIRF0IL  SPANWISE  COORDINATE  »*C13,5 
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1 * SFCTION  LIFT  COEFFICIENT  R*2E13,5) 

«1U  FO«mat  mh  ,/,1m  , /, 3x*AI RFOI L PRESSURE  COEFFICIENTS,  UPPER  ■ *) 

« 1 5 FORMAT  (3X1  OF  1 3,5) 

*tfc  FORMAT  MM  , / , 1 M , /, Sx* AIRFOIL  PRESSURE  COEFFICIENTS,  LOWER  ■ •) 

*17  FORMAT  MM  , / , 1 M ,/,*  UNSTCA0V  FORCE  C0EFFICIFNT9  (PER  UNIT  PITCH 

1  ANGLE  IN  RADIANS)*) 

*18  FORMAT  MH  ,/,3x«LIFT  p*2E  13,5 , / , 3X*MOMENT  ABOUT  (X«XP)  ■ *2EH,9) 
*19  FORMAT  (IM  , / , 1 h ,/,*  PRESSURE  COEFFICIENTS  h*  »IRFOIl  (PER  UN 

1 1 T PITCH  ANGLE  IN  RADIANS)*) 

*20  FORMAT  Mh  ,/,lH  ,/,*  UNSTEADY  FORCE  COEFFICIENTS*) 

*21  FORMAT  (1H  , / , 1 H ,/,*  PRESSURE  COEFFICIENTS  ON  THE  AIRFOIL*) 

*22  FORMAT  MH  ,/,lH  ,/,*  UNSTEADY  FORCE  COEFFICIENTS  (PER  UNIT  PLUNGE 

1 DISPLACEMENT  NORMALIZED  TO  CHORD)*) 

•21  Format  MH  , /, 1 h ,/,«  PRESSURT  COEFFTCIENTS  On  thF  ATRFOIL  (PER  UN 
1IT  PLUNGE  DISPLACEMFtJT  NORMAL  JZEO  TO  CHORD)*) 

RETURN 

END 

SUBROUTINE  GAMFun 

COMPLEX  PHIuB,GamtE!,GAmtF.GAmFF,FPU*FPL 
REAL  *<eAP,MB 

COMMON  /DELTA/  OY(4JO),OY(ttO),DZ(«0),Axl(«0),AX2(<IO),Bxl(ao), 

1 Bx2(R0),CX(a0),AYl((i0),AY2(<l0),AZH2D),AZ2(20),X(«0),Y(«0), 

2 Z(20),FPU(800),FPL(800),PHluB(®00)#IH,IMl,JM,JMl,t(H,XMl,jw, 

3 JwPl.JHHl,ITE,ILE,KSPAN,DYBUl,DYBU2#DYBLl#DYBL2,SMALLK»OHCG, 

a NDOUB, CPCPR, TITLE (8), MB, DEL, ALPHA, !TYPZ,IOPT.XP,N«PRT,KPRT 120), 

5 ZSPAN.XCAP.RPAR 

COMMON  /gamma/  GAMTF1  (20),GAMTE(?0),PGFF,r.AMFK20) 

DO  10  I ■ l , * SP AN 

GAMFF(I)«GAMTE) (I)*PGFF*(GAMTr(I)-GAMTCl(I)) 

CAMTE1 (I)RGAMTI (I) 
to  continue 

RETURN 

End 

subroutine  INITAL 
COMPLEX  PHIIJ8,FPU,FPL 
REAL  *CAP,M« 

COMMON  /DELTA/  DX(a0),DY(«0),DZ(<l0),AXl(«0),Ax2(R0),Bxl(a0), 

1 BX2(R0),CX(«0),AY1 (U0), AY2((I0),AZ1(20),A72(20),X(«0),Y(«0), 

2 Z(20),FPU(800),FPL(800),PHIUB(800),IH,IMI, JH,JM1,KM,XMI,JW, 

3 JWPJ , Jwmi,itE,ILE,RSPAN,DY8U1,DYBU2,DYBL1,DYI»L?,ShALL«,OHEG, 

« NDOUB,CPCPB,TITLE(8),M8,OEL, ALPHA , IT YPE , IOPT , XP, NKPRT ,KPRT (*0 ) , 

9 Z8PAN,KCAP,RPAR 
e CALCULATE  DX,DY  AND  DZ 
DO  10  i«l » TM1 
Dxm«x(iPi)-x(n 
10  CONTINUE 

DO  20  1*1 , JH 1 
0Y(I)«Y(I*l).r(l) 

20  CONTINUE 

00  30  IH,KH1 

oz(n«z(i4i)-z(n 

30  CONTINUE 

C FINITE  DIFFERENCE  COEFFICIENTS 
00  <10  Ip2,IH1 

axi  (n«Dxci.t)/(nx(i>*(Dx(i-n*ox(!>n 

-96- 


CX(I)«,5/D*(i) 

SO  CONTINUE 

cx(n«#5/ox(1j 

on  50  I«2, JMJ 

00  60  I«2,kmi 

to  Jul*'  z<I_1,*f0I,IJ*Dzfi-> 5 J> 

OV^tirlyc  jtJfJ^ovHLl  2,*°Vf  JWM1  > 

8 T AIRFOU.  BOUNDARY  CONDITION 

noPT.FtoArnoPn  TI°* 

00  TO  KaJ ( KSPAN 
00  60  IalLE.lTt 
L«(I»|)*HMtK 

***,,*  Nt“  rUNCT!o«»LlOCPENDENCt"c»»"M0fT*,,'*Ll-'<*t*(I>***)) 

^^T^n^sss,^  «* »*«.. 

•0  CONTINUE 
TO  CONTINUE 
peturn 
End 

SUBROUTINE  PRINT  f NI TERc > 

I »»*(«(,,, B*I(«0) 

l ^f*O)#FPU(«0O),FPL(8O0)  , PHI  IJB  ran  a i r | 0) , AZ2C20) , XC40) , v(<iO) 

comIISJ  /g?5Iv 

R8PAni«K9PAN«|  CZ0) *CAMTEf20),PGFF,6AHFFC20) 

??r}°  5*1  ^a^AN.KSPAiNl 

iw&s 

IJK.MRMTERJMTJ*  WTBRI, •!,),) 

PHI(IJK)«PHI(IJK).p4WT 

L«MP6(ILE.2I.JM4JW 

IP«(ILE»2)*KM*K 

PMIU8 (LP) «PHI (L) 

UP«ITE*KM*r 
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f 


c 


PHIuB(LP)aPHI (IJK)P2.*PART 
COMPUTE  CP  LOWER  (8)  AND  CP  UPPER  (R) 

00  20  I«ILE,ITE 

LiMP*(i-n*jM*j* 

LP«(I«n*KM*K 

B(I)«-2,*(AX1  ( I ) • (PHI  (LpJM)-PHI(L))4AX2(I)*(PHI(L)-phI(L«JM)) ) 

D(I)a«2a*(Axl(I)*(PHluB(LPtKM)«PHlijB(LP))+AX2(I)*(PHIUB(LP)s 

1  PHIUBCLP-KM))) 

20  CONTINUE 


PHI (IJK)«PMI (IJK) tPART 
WRITE  (6,901)  NIT£RG,K 
WRITE  (6,902)  (0(1),  Iall.E,  I T E ) 

WRITE  (6,903)  NITERG.K 
WRITE  (6,902)  (8(1), !■ IlE , ITE) 

10  CONTINUE 

901  PORMAT  ( 1 H ,/,*  at  ITERATI0N#I5*  AND  K «*I3*  SCALED  PRESSURE  COEFP 
1ICIENT,  UPPER  (ILE  TO  ITE)  ■*) 

901  FORMAT  (1H  ,/,*  AT  IT£RATI0N*I5*  AND  K a*!!*  SCALED  PRCSSURE  COEFF 
1ICIENT,  LOWER  (ILE  TO  ITE)  ■*) 

902  FORMAT  (10E13.5) 

RETURN 

End 

SUBROUTINE  TRI  (I,K) 

complex  b,o,pmiub,pmi,fpu.fpl,p 
real  rcap,ms 

common  /DELTA/  DX(40).DV(40),DZ(40),AX1(40),AX2(40).BX1(40). 

1 BX2(a0),CX(«0),AVl  (40),AY2(40),AZU20),AZ2(20),X(40),V(40). 

2 Z(20) ,PPU(B00),FPL(B00),PHIUB(B00),IM,IM1, JM, JM1,km,kmi , Jw, 

3 JWP1,JWM1, ITE,ILE,KSPAN,0VBUl,DVBU2,DVBLl,DVBL2,8MALLK.OMEG, 

« NDOUB,CPCPR. title (»), MS, OEL, ALPHA, ITVPE,IOPT,XP,NKRRT,KPRT(20), 

S ZSPAN.KCAP.RP4R 

COMMON  /COEFF/  A(40),B(40),C(40),D(40),PHI(11500) 

MPbIM*JM*(K«1) 

DO  10  KKO.JMl 

jRJMl«KKtJ 

PRA(J-1)/B(J) 

8(J-1)»B(J-1).P*C(J) 

D(J*1)bD(J«1)«P*D(J) 

10  CONTINUE 

H«MPt(l.l )*JM 
PHI(M»?)«0(2)/B(2) 

DO  20  J«3,JM1 
LRH»J 

PHI(L)»(0(J)-PHI(L-1)*C(J))/B(J) 

20  CONTINUE 

RETURN 

END 


SUBROUTINE  WAKE  (U, SMALL*, RH,WAKEIN) 

COMPLEX  PARTI, PART2.PART3,PART9,CKRAU.CKRH,WAKEIN 
REAL  KRH 
DIMENSION  8(12) 

DATA  c /,3T2/ 

DATA  (S(I),I«1,12)  /!,,«, 24 1861 98,2 ,7918027, •29,991079, 11 1,591 96, 

1 *271, 43549. 305, 75288, 41. 18363,. 545, 98537, 644, 78155,-328. 72755, 

2 64,279511/ 
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CALCULATE  II  F 09  WAKE  INTEGRAL 
IF  (SMALLK.EO.O.)  GO  TO  15 
PART2«CMPL*(0.#0.) 

PART3»PART2 

PAPT4>PART2 

AU»ABS(U) 

SUMQRTU  .*11**2) 
krh«smallk*rh 

CKRH«CMPLX(0,#KRH) 

EKRAU«CEXP(CMPLXfO, ,-KRH*AU)) 

papti«-au*fkrau/su 

DO  10  !>1 « 12 
AMbFLOAT ( I»1 ) 

PAPT2»PART24B(I)*EKRAU*EXP(»C*AM*AU)/(C*AMACKRM) 
10  CONTINUE 

IF  (U.GT.O.)  GO  TO  30 

PART3«*PARTl 

PARTttaPARTP 

DO  20  III. 12 

AMBFLOAT(I-l) 

PART«»PARTR-9(I)/(C*AM*CKRH) 

20  CONTINUE 
JO  CONTINUE 

PART?aPART?*CKRH 

part«»-part«*ckrh 

WAKElN«PARTl*PART2*2,*RrAL(PART3*PARTa) 

RETURN 
15  CONTINUE 

WAKEIN»CHPLX(l,-U/SQRTCl .fU**?)#©,) 

RETURN 

END 
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